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ABSTRACT 


There are numerous numerical methods for solving different types of partial dif- 
ferential equations (PDEs) that describe the physical dynamics of the world. For 
instance, PDEs are used to understand fluid flow for aerodynamics, wave dynam- 
ics for seismic exploration, and orbital mechanics. The goal of these numerical 
methods is to approximate the solution to a continuous PDE with an accurate dis- 
crete representation. The focus of this thesis is to explore a new Discontinuous 
Galerkin (DG) method for approximating the second order wave equation in com- 
plex geometries with curved elements. We begin by briefly highlighting some of 
the numerical methods used to solve PDEs and discuss the necessary concepts to 
understand DG methods. These concepts are used to develop a one- and two- 
dimensional DG method with an upwind flux, boundary conditions, and curved 
elements. We demonstrate convergence numerically and prove discrete stability of 


the method through an energy analysis. 
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CHAPTER 1: 
An Introduction to Solving Partial Differential 
Equations 





There are numerous numerical methods for solving different types of partial dif- 
ferential equations (PDEs) that describe the physical dynamics of the world. The 
goal of these numerical methods is to approximate the solution to the continuous 
PDE with a discrete representation [1]. Three notable methods are Finite Difference 
Methods (FDMs), Finite Volume Methods (FVMs), and Finite Element Methods 
(FEMs). Two common FEMs are Continuous Galerkin (CG) and Discontinuous 
Galerkin (DG), each of which comprises an element-based approach to solving a 
set of equations. The main focus of this thesis is to explore a new DG method, in- 
troduced by Appel6é and Hagstrom [2], for approximating the second-order wave 
equation. However, my research differs from that of Appel6 and Hagstrom by 
using a nodal form of DG with provable stability on curved elements with complex 


geometries. 


1.1 Finite Difference Method 


Finite difference methods are one of the simplest and oldest methods for solving 
partial differential equations [1]. Furthermore, it is arguably one of the most used 
methods for discretizing partial differential equations [3]. There is an enormous 
amount of published information about FDMs in various scholarly journals and 
books. This section describes only the basics of FDMs, and the interested reader 
is directed to, for instance, Gustafsson, Oliger, and Kreiss [4]. Finite difference 
methods focus on approximating the derivatives of the solution directly at a set of 
points in a domain. In terms of the calculus of finite differences, we are looking to 
approximate the derivatives by linear combinations of the function values along a 
grid of points [5]. 


One method of constructing the discretization is accomplished by a Taylor series 


expansion on a selected set of equally spaced points (i.e.,...< Xj-1 < Xj < Xi41 <...). 


For instance, suppose you wish to approximate the derivative of a function evalu- 
ated at x; using grid points x;_1 and x41. By conducting a Taylor series expansion 
of the neighboring values around x;, we can construct a linear combination to get 


an accurate approximation of the derivative. 


For example, let us consider the one-dimensional advection equation, 


Ou ou 
—+—=0, 1.1 
ot ox ey) 
where u = u(x,t) is the solution. We assume an appropriate set of initial conditions 
and boundary conditions for u(x,t). Using the center stencil, we can Taylor expand 
in space around x;, with grid spacing Ax, to find a derivative approximation for gu 


where u; = u(x;,t) [3]. Starting with 








Win = Uj + AX + a +O(Ax), (1.2) 
Ou; Ax? d2u; 
Uj_1 = Uj—-AX e i a + O(Ax’), (1.3) 


(uj4, and uj_1 are not boundary points) and taking the linear combination of (1.2) 
and (1.3) yields 


OU; Ui41—Uj-1 2 

= = ——— + OA"), 1.4 

ox 2Ax ae 
which can be substituted into (1.1) to yield a system of ordinary differential equa- 
tions. Lastly, we can solve these ODEs computationally through a numerical 
integration scheme, such as Runge-Kutta, where 


Ou; Uj Win —Ui-1 
See ee 1. 

at Ox Dax ee) 
Another method for constructing difference formulas is to build a polynomial 
interpolant through the given grid of points, such as using Lagrange polynomials, 
and evaluating the derivative of this polynomial. For instance, using the grid of 


points [x;_1,%;,Xj+41], we have the interpolant 


1 


Po(x) = ye UiseL(X), (1.6) 


k=-1 


where L;(x) is the Lagrange polynomial (see Section 2.1). The derivative approxi- 


mation is then 


1 
Ou; Uist —Ui-1 
Diet i Ly (Xi )Ui+k © — BiRge —o (1.7) 
which is substituted back into (1.1). The Taylor series approach and the polynomial 
approach are equivalent when the maximum number of error terms are eliminated 


from the linear combination. 


1.2 Finite Volume Method 


Finite volume methods differ from FDMs in that instead of seeking pointwise 
solution values, we are seeking the average values over the elements using approx- 
imations. Depending on the dimension of the problem, examples of elements are 
cubes, triangles, or intervals, all organized in an unstructured fashion across the 
physical domain [1]. From the average values, the solution is reconstructed in or- 
der to evaluate a numerical flux to tie neighboring elements together and produce 
an approximation for the entire system. There are numerous FVMs and here we 
only discuss the very basics; interested readers are directed to LeVeque [6] for more 


information. 


Let us again consider (1.1) on a domain that is represented by a collection of 
elements. For this example, let us use elements on a one-dimensional domain and 
each element is denoted by the index 7. Let us further define a different set of 
notation with ou =u; and gu = uy. Therefore, (1.1) becomes u; = —uy. Generally, the 


average value of u on an element is 


Rj 
Q,(t) = =f u(x,t) dx (1.8) 


where Ax, in this case, is the distance between the left and right boundaries of 
the one-dimensional interval and Q;,(t) is the spatial average value of u on the ‘he 
element. If we were using cube-shaped elements, then Ax would be the area of the 


cubic element. If we take the derivative of (1.8) with respect to time, we find 








dQ; 1 x 
= — d 1.9 
dt Ax is ie eo 
and since u; = —u,, we can substitute —u, to obtain 
dQ; i faa 
= —-— d 1.10 
dt Ax f ree am) 


and integrate (1.10) to find an equation to build an approximation to the system: 


dQ; 
dt 





- ~~ RN oUEL OI: (1.11) 


Equation (1.11) is an update equation for the average value of the solution in 
the i' element, where the time derivative of the average value changes through 
the fluxes of the left and right boundaries of the element. There are a variety of 
numerical flux methods to accomplish this approximation. One method is the first 
order upwinding method, meaning that we select the average value of the interval 
always on the “upwind” side. This topic is discussed in further detail in Chapter 3. 
Ultimately, the FVM uses the average element values from across the domain to 
reconstruct an approximation for the system. 


1.3 Finite Element Methods 

Just like FDMs and FVMs, there are many different FEMs used to solve PDEs. 
Similar to the FVM, the physical domain is mapped into a reference domain of 
various sized geometrically flexible elements known as the grid. After the grid is 
built, one of the many different FEMs is applied to solve the system. As discussed 
at the beginning of this chapter, the FEM that is the focus of this paper is known as 
a Galerkin method. There are two main categories of Galerkin methods, Bubnov— 
Galerkin and Petrov—Galerkin [7]. With Petrov—Galerkin, the test functions and 


basis functions are different, where with Bubnov—Galerkin, the test functions and 
basis functions are the same and reside in the same space [1, 7]. We discuss 
basis functions and test functions further in later chapters. For this paper, we be 
focusing on Bubnov-Galerkin methods, often called simply Galerkin methods [7]. 
As discussed before, the two main categories of Galerkin methods are Continuous 
Galerkin (CG) and Discontinuous Galerkin (DG). Both methods use an elemental 
approach to solving the system in integral form, but the difference between DG 
and CG is mainly whether the continuity between the elements is enforced strongly 
or weakly. In the following two paragraphs, let us lightly touch on both methods 
to give the reader a small insight into each method before focusing the rest of the 


paper entirely on element-based DG. 


When using CG, we build a computational domain subdivided into various-sized 
elements, similar to FVM. However, within each element we use specially selected 
degrees of freedom. For example, in the one-dimensional sense, the degrees of 
freedom could be a grid of points across the cell. There are many different ways 
to build this grid of points, such as using Legendre—Gauss points, equally spaced 
points, as discussed in FDM, or some other combination or method. For the purpose 
of this research, we use Legendre—Gauss—Lobatto (LGL) points. These are discussed 
in further detail in later chapters, but they are a set of points that cluster toward 
the boundaries of each element and include points at the boundaries. Discretizing 
the integral form of this equation yields a semi-discrete scheme with elemental 
test functions [1]. These test functions are continuous across the domain in CG 
since we enforce the element boundary points to be equal between neighboring 
elements [3]. Each element has its own set of LGL points; however, using CG 
enforces the solution at the element boundaries to be continuous and this yields a 
global mass matrix. This global mass matrix is for the entire domain and must be 
inverted to solve time-dependent problems; this can be computationally expensive, 


depending on the type of problem [1]. 


Generally, DG is similar to CG in requiring numerical interpolation and integration 
as well as building the same computational domain; however, one of the main 


differences comes from the elemental boundaries. In CG, we require the boundaries 


to be continuous, while in DG they are discontinuous and the continuity is enforced 
weakly. The domain in DG is still represented by a collection of elements; the 
union of these elements is accomplished through a numerical flux similar to FVMs. 
Discontinuous Galerkin still uses the same space of basis and test functions, similar 
to FEMs, and each element boundary has its own set of degrees of freedom [1]. 
Therefore, the solution is typically represented by a set of piecewise polynomials 
that are discontinuous at the element boundaries. The numerical flux, which is 
discussed in greater detail later, resolves this discontinuity to assist in finding 
the final solution. Furthermore, the mass matrix is constructed locally instead of 
globally and this allows it to be inverted at a reduced computational cost, yielding 
a semidiscrete scheme that is explicit [1]. Discontinuous Galkerin is the method 
used for the following study. 





CHAPTER 2: 
Discontinuous Galerkin Foundation 





In this chapter, the information and terminology maybe confusing to the reader 
who has no knowledge of Galerkin methods. However, just like pulling out a 
road map in a foreign country, the intent is to show a DG road map and apply 
these concepts to a one-dimensional and a two-dimensional problem in subsequent 
chapters. Initially, this road map may speak a different language to the reader, as 
expected for anyone new to this numerical method; however, by the end of this 
chapter the intent is for the DG foundation to be clear and understandable. 


2.1 Interpolation 

In order to fully implement DG, we first need to construct the building blocks. 
This leads to polynomial interpolation [3]. Nodal interpolation is the construction 
of an Nt! degree polynomial, f(x), that is equal to a given function, f(x), when 
evaluated at a set of x; points, with i =0,...,N, that is fiy(x;) = f(x;) [7]. There are 
many different interpolation methods and our focus for this section is Lagrange 


interpolation. 


Let us explain Lagrange interpolation through an example. We are going to ap- 


proximate the following Gaussian function in one dimension, 
f@=e™, xe[-1,41]. (2.1) 


We begin by generating a grid of LGL points across the given domain in one 
dimension. As discussed in Section 1.3, LGL points are a specially selected set of 
points that includ the boundaries and cluster toward the ends of the domain. This 
clustering of points toward the ends of the domain assists in avoiding the Runge’s 
phenomenon during the approximation, which is the oscillation of an interpolation 
near the boundaries. This phenomenon is evident when equally spaced grid points 
are used with high polynomial orders. Therefore, clustering points toward the ends 
of the domain helps improve the interpolation [8]. After the grid of LGL points is 


generated, we construct the Lagrange basis functions, L;(x), by the equation 


L(x) = ype £36.38 (2.2) 
1 L J (xj —x)) y PeeepoNy ¢ 
j#i 





where x; are the LGL points defined on the domain [—1,+1]. Using (2.2) for all 
the LGL points in the domain generates a set of N“” order polynomial basis func- 
tions associated with each point x; in the domain. With these basis functions, the 


approximation fy(x) is 


N 
fue) = V° Li) fi (2.3) 
i=0 
where fj = f(x;) fori =0,...,N. Moreover, since (2.2) has the cardinal property 
1 fori=j 
Ley . 2.4 
i(e}) { fori # j = 


we have the interpolation property f; = fy(x;). In the end, fyy(x) is an N' order 


polynomial approximation for f(x). 


We now return to finding an approximation for representing (2.1). Figure 2.1 
shows the polynomial interpolation of (2.1) using N = 2, 4, 8, and 16 with N+1 


LGL interpolation points. In Figure 2.2, we consider the error norms 


50 in a 
llelln = a, (2.5) 
ae | F(x) | 
eae — f(&))2 
Ye Gepe 
max,<r<50 | f(x) — f (Xp) | 
max1<x<50 | f (Xx) | 


lle ll = (2.6) 


II € [loo (2.7) 


where we have used an equally spaced grid of 50 points, %,, within the domain 
[-1,+1]. Figure 2.2 displays the error norm of N = 1 through N = 40 polynomial 


interpolations. 


As you can see in Figure 2.1, the 16""- order polynomial interpolation does a good 
job approximating the function as compared to the exact solution (denoted by * to 
differentiate from the various interpolations). In Figure 2.2, we show the error in 
the interpolation vs. the N‘" order polynomial interpolant based off of (2.5), (2.6), 
and (2.7). As the figures show, the approximation gets to machine precision around 
the 32” order polynomial. Lastly, the reason that Figure 2.2 has various plateaus 


within the convergence rate is because f(x) is an even function. 
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Figure 2.1: Lagrange Interpolation of e* at Figure 2.2: Approximation Error based 
different sets of LGL points for various poly- off (2.5), (2.6), and (2.7) measured at poly- 
nomial orders of N. nomial orders N = 1 through N = 40. 


2.2 Integration 

The next building block for DG is the ability to conduct numerical integration. 
When performing numerical integration, also known as quadrature, the analysis is 
similar to that of interpolation. Integrating the Lagrange interpolation (2.3) gives 


the integral approximation 


+1 +1 _N N +1 
{ fiv(xdx = { 2 Lied = ) f(x) { L(x)dx. (2.8) 


Defining the quadrature weights 


+1 
wee { Lj(x)dx (2.9) 


1 


gives the numerical integration method [8] 


+1 N 
ax = iJi- Al 
‘a fu(x)dx Dwih (2.10) 


These quadrature weights only have to be precomputed once with the LGL points 
and not at each evaluation of a particular integral [8]. Figure 2.3 is the comparison 





Error 
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foo} 
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to! + 


10°12 L 














Figure 2.3: Numerical Integration Error evaluated with the Euclidean Norm vs. the N‘ order 
quadrature. 


: : : ‘ ; 1 
of the numerical integration approximation to the actual value of f ¥ eA) dy x 
0.882081. The Euclidean Norm is used for the calculation of the error due to the 
integration of f(x) producing a scalar. As you can see, the error decreases to 


machine precision around the 18"" order quadrature. 


10 


In comparing Figure 2.2 to Figure 2.3, we see that when using the same order 
polynomial, integration is more accurate than interpolation. In general, using 
N+1 points, we can construct an N“" order polynomial. Since we fixed the LGL 
points at the boundaries, -1 and +1, there are only N—1 LGL points left to choose. 
We also have N +1 quadrature weights to choose. Therefore, we have N — 1 points 
and N +1 weights, which means there are 2N degrees of freedom. We can thus set 
the degrees of freedom so that we integrate a polynomial of order 2N —1 exactly 
(as this polynomial has 2N coefficients) [3]. For example, we can exactly integrate 
a fifth degree polynomial using N + 1=4 LGL points; even though we can only 


construct a third-order interpolating polynomial using these points. 


2.3 Concept of the Mass Matrix 

After discussing interpolation and integration, let us now consider the concept 
of the mass matrix. The mass matrix, denoted by M, is used for integration. 
We discuss the integration of the product of two interpolants exactly and inexactly. 
Using exact integration results in a full mass matrix, while using inexact integration 
results in a diagonal mass matrix. For example, suppose we wanted to integrate 
f(x) with a test function, (x), 


+1 
{ ; f(x) (x)dx, (2.11) 


over the domain of a single element from [-1,+1]. Test functions are defined in 
further detail in Section 3.3, but for this section all we need to consider is that (x) 


is a function that resides in the same space as f(x). 


Suppose both (x) and f(x) can be numerically approximated with Lagrange poly- 


nomials. Then 


N 

W(x) * pre) =) POx)Li(x) and (2.12) 
. 

Fx) = fulx) =)" faedLin). (2.13) 
i=0 


11 


This leads to the integral approximation 


+1 +1 +1 N N 
[, foopenae~ [ pucopnonde= [YY foatLacateoiverpar, 


-1 ‘=0 j=0 


— 


which is further simplified to 


N N “4 
YY ve ie Li(x)L jonas | fen) = TM, (2.14) 


i=0 j=0 


~~ 


mm. 


where 
+1 
{ L,(x)L j(x)dx = Mj; (2.15) 
=f 


and Mf; is a full matrix per element that can integrate an order 2N function exactly. 
Furthermore, y and f are the evaluation of (x) and f(x) at the LGL points stored 


in column vectors: 


Ff (x0) (x0) 
fz fe) 2 eae ve) 


f (xn) W(xn) 


Inexact integration is similar to exact integration. Instead of numerically approxi- 
mating the functions separately with Lagrange polynomials as seen above, we first 
multiply the function and test function together at the LGL points and produce a 


new interpolant for numerical integration: 


+1 _N 


im feb (x)dx i 2 fedveay (x)dx = Ys (2.16) 


Once again L(x), under integration, generates the weights, wj;, as seen in (2.9), 


and w; is only computed once at the LGL points. These weights can be stored 
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in a diagonal mass matrix, called M, and used in an one-dimensional problem. 


Therefore, (2.11) can now be written in the matrix form below: 


+1 
{ _ Owe de yy Mf. (2.17) 


As discussed in Section 2.2, using this method with LGL points allows us to inte- 
grate a polynomial order of 2N —1 exactly. Throughout the remainder of this thesis, 
we only consider inexact integration. This is because the diagonal structure of the 
mass matrix greatly simplifies the implementation of the method. Additionally, the 
Jacobians and surface Jacobians are easier to handle when moving into multiple 
dimensions. We discuss more about Jacobians and surface Jacobians in subsequent 


sections and chapters. 


2.4 Concept of the Differentiation Matrix 

The differentiation matrix has a purpose similar to that of the mass matrix, but it is 
obviously used for differentiation. Like the mass matrix, the differentiation matrix 
is found using Lagrange polynomials by 











aL j(x;) 
ij = va : (2.18) 
d fn (xi) ae 
e 4 i= Lat (2.19) 


where it should be evident from (2.19) that Df approximates af at all the nodes. 
As you can see, D is simply the derivative of the Lagrange polynomials evaluated 
at the LGL nodes. You should further notice that D is a full matrix. 


Let us examine how the differentiation matrix is used for discretization. Consider 
the integral 


+1q oe 





W(x)dx, (2.20) 
where 1(x) is a test function that resides in the same space as f(x). Once again, the 
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functions are numerically approximated with a set of Lagrange polynomials 





N 

d(x) © Y pexLilx), (2.21) 
= 0 

oe x yf Eo i(X) (2.22) 


leading to 


aoe +1_N N 
[. Pocars [7 YY soa| ——e 09] perp, (2.23) 


i=0 j=0 


which further simplifies to 





N WN 
VL ve | (ee aE Or, io fc)) » "MDF. (2.24) 


As you can see, (2.18) is nestled right in the middle of (2.24). Once again, both p 
and f are the evaluation of (x) and f(x) at the LGL points. 


2.5 Change of Variables 


Mapping from a physical space to a computational reference space requires a 
change of variables. In this case, we are going to map from xX € [Xy,Xy41] to € € 
[-1,+1], where n =0,...,N. Furthermore, a Jacobian arises when conducting this 
change of variables [9]. In one dimension, the Jacobian is simply A where the 
element size is h = Ax for an equally spaced grid of elements. This one-dimensional 
Jacobian is used extensively in Chapter 3 for discretization and is annotated by 
J= ft In two dimensions, it is a little more complicated and is discussed in greater 
detail in Chapter 4. However, the concept remains the same for one dimension and 


two dimensions. 


In one dimension, we want the computational reference domain to be from € € 


[-1,+1] for two points from a linear element, x; and x,+1. This concept is depicted 
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in Figure 2.4. 


= 
Xn antl é9=-1 €,=+1 


Figure 2.4: Physical Space to Computational Space 


The change of variables for Figure 2.4 is obtained by a linear combination 


x(€) = (SS), + (* : 2 sna (2.25) 





that allows us to approximate the coordinates of the element. For example, if € =—1 
in (2.25), then that maps to x, and if € = +1, that maps to x,41. Taking the derivative 
of (2.25) yields 


dx = (Xn+1—Xn) _ h 


77 : 5 (2.26) 


and dx = Adé [3]. Therefore, using (2.25) and (2.26) for the change of variables for 


integration gives 


X41 +1 h 
fo fade =f Fptotende, (2.27) 
where xX, and Xx,+; are the left and right boundaries for the physical element. 
Constructing this change of variables for each element is a key foundation for local 
element based Galerkin methods. Instead of building different matrices for each 
element and solving them individually, we can now construct the required matrices 
for a reference element and then use the metrics terms to scale the reference element 
to the physical space [3]. Metric terms are discussed in greater detail in Chapter 
4. As for using the differentiation matrices, D, the change of variables requires the 
chain rule, which creates a ae = c multiplied by a = A cancelling these terms out. 


The equation below is a visual example for the above sentence showing the change 
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of variables for (2.20): 


“yon fOace "ey ESIC te 0 ‘LOO ae = yoy 


Xn 


This section assists in the construction of the one-dimensional grid further investi- 


gated in Section 3.1. 


16 





CHAPTER 3: 
One-Dimensional Discontinuous Galerkin 





Let us start working through the Discontinuous Galerkin method for the linear 


acoustic wave equation in one dimension. Consider the equation 


Op 2x2 d 

= = Vp, xEOCR’, t>0, (3.1) 
ot? 

where p is the pressure, c > 0 is the wave speed, and d is the spatial dimension. On 

the boundary, represented by 0Q, we impose periodic boundary conditions. We 

are now going to split this second-order wave equation into two equations with 


the auxiliary variable, w, equalling the time derivative for pressure [2]: 


Ow mee) 
es =¢ Vv P, (3.2) 
op 


3.1 One-Dimensional Grid 

Figure 3.1 shows an example of a one-dimensional grid of equally spaced elements, 
with polynomial order N = 6 and three elements, covering © = [—1,+1]. The LGL 
points, depicted by red x symbols, are used for the interpolation points per element. 
Polynomial order, N, is defined as the maximum order of a polynomial that can be 
represented exactly on each element. Obviously, the grid varies depending on N 
and number of elements (Ne). As you can see, the LGL points are not evenly spaced 
across the individual elements and cluster towards the boundaries of each element. 
There are N +1 LGL points per element, where each element is represented by Q), 
with j = 1,...,Ne. Remember, this is a one-dimensional grid with Q = [-1,+1] 
and each element (Q/;) is mapped to é € [—1,+1] through a change of variables, as 
discussed in Section 2.5. This one-dimensional grid is simple; however, in two 
dimensions, the grid is a little more complicated and is discussed in Chapter 4. 
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-0.333 +0.333 


Figure 3.1: One-Dimensional Grid with N = 6 and Ne =3. 


3.2 Variational Form 

Understanding (3.2) and (3.3) motivates the variational form that is needed for us 
to build the set of equations for DG. This section explains the process of finding 
the variational form of (3.2) and (3.3). Furthermore, the discretization of these 
equations yields the necessary equations for the DG numerical approximation. 


First, we need to focus on finding a variational form. 


3.2.1 First Variational Equation 


Equation (3.4) is the most trivial of the three variational form equations. If ep =o, 


it follows that 
Op | 
—-w]=0. (3.4) 
Lig 


Using (3.4) in Section 3.4 assists in finding the needed equations for the DG ap- 


proximation. 


3.2.2 Second Variational Equation 
Given that a0 = c°V*p, by multiplying this equation by a test function, , and 


F) 
{, | y (= - vp) -0. (3.5) 


We define test functions in further detail in the following section. By conducting 


integrating over Q; yields 


integration by parts on (3.5) and inserting a numerical flux at the boundary of the 
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element, we are able to find the second equation in the variational form to be 


i (y5¢+2vy- ve)=2 3 wn-(Vp)’. (3.6) 
OF 


For notation purposes, a * term in the equation identifies the numerical flux. Using a 
numerical flux is a numerical technique for coupling the elements that is commonly 
used in DG. For example, (Vp)* is the coupled numerical flux computation for the 
gradient of the pressure at the boundaries (0Q/;) for neighboring elements within 
the domain. An explicit computation of the numerical flux is discussed in greater 
detail in Section 3.5. 


3.2.3 Third Variational Equation 
In order to find the last variational equation, we are going to multiply the gradient 


of (z - «| with the gradient of a test function (Vw): 


op = op _ 
[vers -eo)= [vers ff veve=0 (3.7) 


Let us focus on the last portion of (3.7). By conducting integration by parts we find 
{ VV = i (Vp-n)w- { (Vp) a. (3.8) 
Q; dQ; Q; 


By introducing a numerical flux to the boundary term of (3.8) and substituting it 
into (3.7), we find that 


d * 
fF [comes [eens 


J 


Let us focus now on the last term from (3.9). Conducting integration by parts again 


on the last portion of (3.9) yields 


je Jo= fo, (Vip-n)w— I, VipVw. (3.10) 
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Next, we substitute (3.10) back into (3.9) and combine like terms to conclude the 


following: 
op : 
VyV—-— (Vp-n)(w* -—w)- VyVw = 0. (3.11) 
Q; at Jaa, Q; 
Finally, from (3.11), we arrive at the final variational form 
op : 
VyV| —-w]= (Vw -n)(w -@). (3.12) 
Q; ot aQ; 


Overall, equations (3.4), (3.6), and (3.12) together are the variational form of (3.1) 


that we discretize using DG. 


3.3 Test Functions 

Using Discontinuous Galerkin, we are looking for p and w from (3.2) and (3.3) 
that satisfies the variational form found in Section 3.2 for all piecewise polynomial 
test functions. Since we are working to build an approximation to a function with 
piecewise polynomials, our space is the space of N''-degree piecewise polynomials, 
with the objective of solving for a piecewise polynomial that represents the solution. 
For example, suppose we wanted to find a constant approximation for f(x), such 
that the error was orthogonal (in an integral sense) to all other constants. That said, 


we want to find a € R such that 

+1 +1 

ih YGF = i _ (Ola Fax = 0, 
for all CE R. If f(x) = x”, then the problem is to find a such that 
+1 1 +1 
i (Ca - Cx?)dx = |cax — 5cr| =0, 

-1 3 = 

which is simplified to 


2Ca-5C=2C(a-=)=0. (3.13) 


Therefore, a = $ is an approximation for x? found using a 0" degree test function 
(C). Expanding this concept to DG is essentially the same, except we are using N!"'- 
degree piecewise polynomials as the test functions to find a N“" degree polynomial 
approximation for the solution on each element. We use this concept for the 
remainder of this paper and Section 3.3.1 further highlights the use of test functions 
for the discretization in Section 3.4 and Section 4.2. 


3.3.1 Test Functions in Discretization 

As discussed in the previous section, we are using N“"-degree piecewise polynomi- 
als as test functions for DG. In the following section and chapter, we discretize the 
variational form to find a discrete approximation for the continuous equation in 
one dimension and two dimensions. These discrete approximations, with respect 


to y, all have the similar form of 
yy" (Ap - Bw) = 0, (3.14) 


where A and B are matrices and p and w are the solution vectors. Since (3.14) holds 
for all y, then Ap— Bw must equal zero. This concept holds true for Section 3.4 
and Section 4.2. Therefore, in what follows, wy is not listed in the final discrete 


approximations. 


3.4 One-Dimensional Discretization 

In this section, we derive discrete versions of Equations (3.4), (3.6), and (3.12). The 
idea of discretization is to replace a continuous equation with a consistent discrete 
approximation. As you have probably noticed, with Galerkin methods we use the 
variational form for discretizing equations [3]. We focus on developing a spatial 
discretization of the equations and time is integrated separately with a Runge-Kutta 
method. We isolate the time derivatives and combine the discrete forms of (3.4) 
and (3.12) into one discrete equation [10]. As discussed in Section 2.5, the change 
of variables requires the Jacobian of J = u for an equally spaced grid of elements. 
As a reminder, the concept from Sections 3.3.1 is applied to the discretization in 


Section 3.4.2 and Section 3.4.3. We discretize each equation individually. 
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3.4.1 First One-Dimensional Discretization Equation 
Using the concepts established in Chapter 2, the discrete version of (3.4) is 


d 
Myst ~1'MJw = 0. (3.15) 


Within (3.15), 1 is a vector of ones and is needed because in (3.4) we are integrating 
against the function (x) = 1 [10]. Additionally, 7 and w are column-vector ap- 


proximations of the solution at the grid points. This is similar to f from Section 2.3. 


3.4.2 Second One-Dimensional Discretization Equation 
The following is the discrete version of (3.6): 


d op\" op\" 
Me +c2J-1DTMDp = 2Jen| | — eof 2) J. (3.16) 
dt JE} 0 
The column vectors ey and eg consist entirely of zeros except for the last “row” 


in ey and the first “row” of ep being ones. Using ey and eg is a way of ensuring 


* 


that the first and last portion of information from the numerical flux of (32) and 
N 


(3) is used for the calculation. The numerical flux is discussed in greater detail 
0 
in Section 3.5, but it is essentially a method for coupling the solution on either side 


of an element boundary. Additionally, discretization of 


[ vo-vp 


J 


from (3.6) requires two differentiation matrices; D! takes the derivative of the test 
function, M is the integral, and D is the derivative of p. This accounts for the D) MD 
portion of (3.16). 


3.4.3 Third One-Dimensional Discretization Equation 
The following is the discrete version of (3.12): 


dp 


J-'D'MD—- ~J"'D™MDo = J"'D' ey (wy — ew) -J-'1D"e9 (we, w). (3.17) 
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Once again, the {, Vv . ve and {, Vv -Vw from (3.12) are accounted for by the 
J J 


D'MD portions and w* is computed through the flux at the element boundaries. 


3.4.4 Combination of One-Dimensional Discretization Equations 
Now that we have three discrete equations, our goal is to get the problem into the 
form of two equations and two vector unknowns. Combining (3.15) and (3.17) 


yields 


(J-'D™MD+ 1M]) 2 =(J-'D™MD + IMJ) w+J"'D'en (wy - en) (3.18) 
—~J-'D eo (cw —ejw). | 


For notational purposes, 11 is a necessary square matrix of ones of size (VN +1,N +1) 
because (3.15) is a scalar. Additionally, letting ( J-‘D'MD+1M]J ) be equal to the 
variable Mj and simplifying (3.18) produces 


dp a * * = 
Mi aes Miw+IJ Ip? (enw, — eowp) —J Uy! (enexw — evey@). (3.19) 
For the second equation, by rearranging (3.16) we can isolate te as seen in the 
following: 
dw 27-1pT 27-1 op\’ op\ 
a D*MD —J] —eo|—| |}. 2 
MJ, cJ MDp+c*J en\5E}, “laE |, (3.20) 


Equations (3.19) and (3.20) are the two systems of ordinary differential equations 


that we solve numerically using a Runge-Kutta method. 


3.5 Numerical Flux 

When using a DG method, we have to account for the discontinuity that exists 
between the elements. Meaning, when each element is represented by a polynomial 
for w and p, sampling w and p at the boundary of neighboring elements, as shown in 
Figure 3.2, yields different results [3]. The numerical flux is a method for enforcing 
(approximately) continuity between the elements. There are multiple flux methods, 


but in what follows we only discuss the central and upwinding fluxes. 
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Figure 3.2: One-Dimensional Element Boundary 


The central flux is the easiest method to understand, but it does not necessarily 
produce the best results. The central flux is the average of the values at each 
boundary element. In Figure 3.2, the neighboring elements, 0; and Qj+1, have 
different approximations at the boundary point (BP). For notation purposes, the 
superscript (+) represents the left element’s gridpoint and the superscript (—) repre- 
sents the right element’s gridpoint. Furthermore, as in Section 3.4, the superscript 
(+) denotes the numerical flux term. This notation comes into further use in the 


following section and chapters. The central flux, for both w and ae is defined by 


the following: 
* 1 - + 
a= 5 +w'), (3.21) 
dp” 1, (op\ (ap i 


Using the central flux is a useful beginning step when coding the flux into any 
algorithm. The central flux is simple and still achieves convergence when coded 
correctly. However, central fluxes often have suboptimal convergence rates and 
can lead to oscillatory approximations due to a lack of dissipation. To improve this, 
we need a numerical flux that does not just average the information at the element 
boundaries, but can account for the physical propagation of information across the 


boundaries. This leads us to the upwinding flux. 


An upwinding flux is a method for computing a flux across the element boundaries 
based on the physical propagation of information. In the acoustic wave equation, 
the waves that are propagating through the physical system can move in both 


directions across the element boundaries. We are looking at the information that 
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is “upwind” of the wave’s direction of motion, hence called an upwind flux. For 
example, for the advection equation, if the wave is moving across the element 
boundary from the left to the right, we select the information from the left element 
and vice versa for the wave moving from right to left across the boundary. However, 
for the wave equation we want to decouple the wave propagation into two one-way 
advection equations to take advantage of this “upwind” concept. 


Since we are working in one dimension, let us take (3.2) and derive the upwind 


flux. Defining the auxiliary variable g = px, we rewrite (3.2) as 
ear pe (3.23) 


Additionally, understanding that py: = pix, we can further deduce that g; = wy and 
write a first order system of equations as 


q|, 1 O}lg & 
By defining U = "I we simplify the system of equations into a matrix vector form: 
q 


0 c 
u,;-AU,=0, A= “or (3.24) 


Solving for the eigenvalues and eigenvectors of A, we find the eigenvalues, A; =c 


and A» = —c, and the eigenvector matrix (V) and its inverse (V~!) to be 


We now want to find a linear combination of the original variables, w and q, whose 


Cc 
V= 
1 





fa and V-!= 
1 





a 
NIF NIF 


solutions are independent of each other. We can do this by diagonalizing (3.24): 


[v-tu] -[v-tav][v-tu]_ =o. 


25 


We define a new set of variables r; and rz, known as the characteristic variables, by 


pov S|" 
12 , 


and let A be the diagonal matrix of eigenvalues 


A, O 
[vtav]=aq| OO. 
0 Az 
We now have a system of first order one-way advection equations that allow us to 
upwind because the propagation of the waves are decoupled across the element 


boundaries by the characteristic variables, r; and rz, as seen with 


aches at Oe (3.25) 
; * P21 0 Ao 2 |. Ol 


By analysis of the function within the element, with f,; and fp being the initial 


conditions for 7; and 12, we find that 
r1(x,t) = fi(xt+Art) = fi(xt+ct), and ro(x,t) = fo(x + Act) = fo(x —ct), 


which further describes the propagation of the wave across the element boundaries 
with 11, flowing right to left, and rz, flowing left to right. Figure 3.3 depicts this 
propagation. We know the direction of propagation and can now find the flux 


values. By letting r; =r, and r5 = ry, we have chosen the numerical fluxes as the 





Pai oer 
QO; L +9- R Oj+1 
rt , 
2. 


Figure 3.3: One-Dimensional Element Boundary Decoupled Wave Propagation 
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upwind values. Knowing that r = V~1U, we can easily find a system of equations 
to solve. Note that 


_ 
ma ~ 
N ae 

[eae ess 

I 

. lee 

YR 

NIF NIFR 


so with if =1, and i = ae we have 
1 1 1 1 
"=—@M@ +-g =— eel er d 3.26 
ge ee ape, oe oe) 
1 1 dys We 
es = cal + ue = ae + 54 3 (3.27) 


We now have two equations and two unknowns where we can solve for w* and q"*: 


Pe ae 1 : 
Pas -e )+5 +9); (3.28) 


ges ; ("+0") + SQ--4"). (3.29) 


Lastly, with q = ae (3.28) and (3.29) are the upwinding flux equations for w and i 
in one dimension. It should be further noted that both central flux equations, (3.21) 
and (3.22), are included in (3.28) and (3.29). The portion of (3.28) and (3.29) that 
include the variable c is considered the upwinding portion and is used separately 
for the upwind flux energy analysis in Chapter 5. 


3.6 One-Dimensional Discontinuous Galerkin Results 
After laying the groundwork, we now apply these concepts to an actual problem 


in one dimension using an upwind flux. We use the exact solution 


p(x,t) = sin(n7x) sin(n7t), (3.30) 


Pew=nn sin(n7x) cos(n7t), (3.31) 


Ot 


on the domain Q = [—1, +1] with periodic boundary conditions where 7 is an integer. 
To build the initial condition, (3.30) and (3.31) are evaluated at t = 0 across a one- 


dimensional grid of LGL points where x € [—1, +1]. These initial conditions are used 


2/- 


within (3.19) and (3.20) and evaluated through a Runge-Kutta five-stage fourth- 
order (RK54) accurate iterative method [11]. Runge-Kutta methods are numerical 
methods that approximate a system of ordinary differential equations. With the 
RK54 method, (3.19) and (3.20) are evaluated through five stages per time-step 
along the grid of LGL points. All of these evaluations are combined to produce a 
fourth order accurate or higher approximation [8]. The following section discusses 
the results of applying a DG method for the above one-dimensional problem. The 


DG implementation can be found in Appendix B. 


The implementation is tested using various polynomial orders and numbers of 


equally spaced elements within the domain listed in Table 3.1. Once again, the 


Table 3.1: Discontinuous Galerkin Tested Information 


Polynomial Orders (N) | Number of Elements (Ne) 
N=4,6,8 Ne = 2,4,8,16 
N=16 Ne =2,4,8 























DG algorithm uses inexact integration for computational speed to ensure an easily 
invertible diagonal mass matrix. Figure 3.4 shows the log of the error vs. log of the 
number of elements for N = 4,6,8 calculated using the global L? error norms. As 
expected for both w and p, Figure 3.4 displays increasing convergence rates as the 
polynomial order increases. Increasing the order of the local approximation gives 
the fastest convergence rates, as apposed to increasing the number of elements, 
due to the fact that the global error is dependent on the polynomial order [1] as 
depicted by 


ll € lz < Chit. (3.32) 


In (3.32), C is a constant that is not dependent on the element size h but does depend 
on the final time, t, of the solution. In this case, N is proportional to g in (3.32). As 


the polynomial order (N) increases, the convergence rates increase. 


Table 3.2 displays the convergence rates for the tested information from Table 3.1 
for w and p. As you can see, as the polynomial order increases, the convergence 


rates increase. Furthermore, Table 3.2 depicts the convergence rates being near or 
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Convergence Plot for w Convergence Plot for p 





5 10 15 5 10 15 
Ne - Number of Elements Ne - Number of Elements 


Figure 3.4: Convergence Rates For N = 4,6,8 at Ne = 2,4,8,16 


better than its associated polynomial order, where the convergence rate for w is of 
order N and the convergence rate for p is of order N +2. 


Most DG methods are expected to be order N, N + 5, or N+1 [1]. However, the 
results yielded N and N +2 when using the same space of functions for w and p. 
In Appelé and Hagstrom’s paper [2] they proved that when p is an N +1 order 
polynomial and w is an N order polynomial you get optimal convergence of N+1 
for this method [2]. In Table 3.2 and 3.3, these convergence rates are an observation; 
a proof would require further analysis. 
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Table 3.2: Convergence Rates for N = 4,6,8 






































Convergence Rates 
N/Ne | Ne=2to4 | Ne=4to8 | Ne=8to16 
N=4(a) 3.2286 4.1370 4.1317 
N=4(p) 3.2103 6.3134 6.8160 
N=6(q) 5.8699 6.0861 6.0456 
N=6/(p) 7.2975 8.9715 7.9343 
N=8 (a) 7.7170 8.0165 8.0291 
N=8(p) 8.9915 10.1413 9.8260 














Table 3.3: Convergence Rate for N = 16 























Convergence Rates 
N/Ne Ne=2to4 | Ne=4to8 
N=16(w) | 14.0437 15,6555 
N=16(p) | 16.5784 17.5193 














When testing higher order polynomials, we can increase the spatial error for the 
given function and decrease the time-step in order to ensure that we can neglect any 
time-stepping errors during testing [1]. Doing this allowed accurate convergence 
rate measurements when testing 16" order polynomials. Figure 3.5 and Table 3.3 
display how even higher order polynomials can be used to approximate the solution 
with a much higher convergence rate. However, using higher order polynomials 
comes with a computational cost with respect to time. Higher orders are much 


more accurate, as seen in Figure 3.5, but the computational time increases. 
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Convergence Plot for w Convergence Plot for p 
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Figure 3.5: Convergence Rates For N = 16 at Ne = 2,4,8 
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CHAPTER 4: 
Two-Dimensional Discontinuous Galerkin 





Moving into two dimensions is obviously more challenging than one dimension. 
The first portion of Chapter 4 focuses on an example two-dimensional grid with 
quadrilaterals and explains the formulation of the metric terms needed for dis- 
cretization. Section 4.2 is the discretization of (3.4), (3.6), and (3.12) in two dimen- 
sions. As discussed in Chapter 2, working with integrals and mapping the physical 
domain into a reference domain requires a change of variables. In two dimensions, 
this change of variables produces a set of Jacobian determinants as well as surface 
Jacobians on each element. Lastly, we apply these concepts to a two-dimensional 


problem on a washer domain with curved elements. 


4.1 Two-Dimensional Grid 

Instead of a one-dimensional grid, which is just a line of elements in Chapter 3, 
moving into two dimensions requires a two-dimensional grid of elements. Building 
the grid requires a change of variables, from the x and y spatial coordinates to the 


€ and 7 reference coordinates [9]. Figure 4.1 is an example of a mapping for 








3 3 
8 
5 

4 2 => 4 x 2 

2 

1 1 

y => it 
x c 


Figure 4.1: Example Mapping of a Quadrilateral Element in Two Dimensions with Degrees of 
Freedom for N = 2. 


a quadrilateral element in two dimensions. When using quadrilaterals for two 
dimensions, we have (N +1)? LGL points in both the rand € directions, as depicted 


uo 


by the red xX symbols in Figure 4.1. The numbered order of the LGL points in 
Figure 4.1 is important for the storage of data. For example, in Section 4.2 the 


solution is stored in column vectors 


W1 Pl 

W2 p2 
@= , P = . ; 

WM PM 


where there are M = (N +1)* LGL points or degrees of freedom. 


For notational purposes, the matrix of Jacobian determinates is annotated by J and 
the matrix of surface Jacobians is annotated by one of the four sides of the element, 
such as Sj; for side one of the element, as labeled in Figure 4.1. Though Figure 4.1 
is an example of elements with straight boundaries, we develop the scheme for 


general curvilinear quadrilateral elements. 


4.1.1 Metric Terms 
As discussed before, we need to conduct a change of variables from (x, y) to (€,7), 


where 


x=x(2,n), y=y(én), (4.1) 


and we assume the inverse mapping exists with 


G = alte Y), if oe n(x, y). (4.2) 


Let us consider a single quadrilateral element, such as the one in Figure 4.1. The 
following explanation for the metric terms is based on Professor F.X. Giraldo’s 
lecture notes [3]. Consider the differentials defined from (4.1) and (4.2). They are 


Oe) OG 08 
ax= BES + pra dé = oy + ay 
_ oy Oy _ on on 
dy = pE4° + ee dn = Ft + ay 
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which can be written in the matrix form 


Ox ox 
dx| \3e 3p a 
= : (4.3) 
‘ E i dn 
HB aL) “ 
dn| |e 3, \l4y 


The Jacobian and the inverse Jacobian, along with their determinants, from (4.3) 
and (4.4) are 


ox ox 
N88 onl ee _ axdy dy dx 
S=la, a | 1s = 5h es (4.5) 
Ei a dEAn JEAN 
0& 0 
2 Aa pial 2 4 déOn Onde 
a eg 
x Oy y y 


Focusing on (4.5), taking the inverse of the Jacobian, J, yields 





if 2 
(D7 = 5) By Hl (4.7) 
OE OE 
which must equal J lin (4.6) and thus 
= s|_1/% - 
x oy|_ =| on on 
E -1% 3] “i 
ox oy 0g 0g 





Equation (4.8) gives the metric relations 


a dy aq __ ay A__ ax aq_ax es 
ox on’ “dx  0& “dy on’ “doy d&’ 
which arise after the chain rule is used to map from the physical domain to the 
reference domain. Lastly, on each element, we store the Jacobian determinants 
and metric terms (4.9) in diagonal matrices, where the storage order is as the LGL 
gridpoint order in Figure 4.1. 


os) 


4.2 Two-Dimensional Discretization 

The beauty of moving into two dimensions with quadrilateral elements is that the 
constructions of the mass and differentiation matrices are relatively simple. The 
evaluation of the surface integral over each element requires the Kronecker product 
of the one-dimensional mass matrix with itself. This allows us to integrate in both 
the € and 7 directions for each element and is annotated by M@M, where M is just 
the one-dimensional mass matrix. Additionally, the differentiation matrices are 


Dz =I@D, 
D, =Del, 


where I is the identity matrix of size (N + 1) x (N +1) and D is the one-dimensional 


differentiation matrix. 


4.2.1 First Two-Dimensional Discretization Equation 
The discretization of (3.4) is very similar to the one-dimensional discretization 


except for the addition of the Kronecker product 
d 
17](M@M)—-1"J(M@M)w=0. (4.10) 


Since we are in two dimensions, 1 is a vector of ones of size (N + 1)” and is needed 


because in (3.4) we are integrating against the function 1)(x) = 1 [10]. 


4.2.2 Second Two-Dimensional Discretization Equation 

The discretization for (3.6) requires much more analysis than (3.4) because of the 
gradients and surface terms. Due to this, we focus on the left-hand and right-hand 
sides separately. Here we explicitly introduce the differential into the integral for 
two dimensions, whereas in Section 3.4, this differential was implicit. Focusing on 
the left-hand side, we see that 


ve . - dw " +1 tl 
{ po +oVp Vp)dA = J(M@M)—- +c ‘i Hi [Vw - Vp] Jdédn, 
Q; -1 -1 
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where Vi) and Vp are equal to the following after the change of variables: 


OE, OE» on, On »| Op 
We ae sulla E By! | On’ 
dE» ap fan, an ay 
ae Fae sallae* E  3y!| an" 


Therefore, Vy - Vp yields eight separate terms: 
LFA ca eLearn clc 
OxdE Axdn||Ax dE AxOdn| |AyIE ODydAn||dyVE dDyoOn 
DE OpPIEIP dE Sp ondpy  Anopd— op dE Op On op (4.11) 


Vy -Vp= 


~ Ox dé ax Oe Ax OE Ax On Axdndx dE Ax dAnOx dn 
dE OpAE Op dE Spandp In opdedop | dE op on op 


* Oy dE Oy dE OydEdy On OydndydE dDydndy On 


Focusing on the first term of the second equation in (4.11), we find 


+1 tl +1 tl 
d&\ dp (d&\ oy _ dé \ Op j2e Ow 4 
i ie (=) 0€ (=) dé TS aa -{ {. Nae 0& Ia Ea aca, Mele) 
where we have multiplied in JJ~! to place a Jacobian determinate with both of the 
metric terms. By substituting in the metric terms from (4.9) and discretizing, we 
produce 


+1 +1 Oy dp dya P) e F) 
J ae “aan = (SPDea) pimem(3 De| 
1 1 1 


-1 ANE An OE (4.13) 


OY: oy 
— pI nt 1 


The same process is completed for each of the eight terms in (4.11) to find the final 


discrete form for the left-hand side of (3.4) to be 


y'J(MeM ae +c" Tp, (4.14) 
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where 


-{ pt|?4 1 29 | 5. pt|?¥ oy 
r=( ot) eM) ;, Dz-D; ml (MM) =, D, 
01) tr (M@M) I: +p er (M@M) “4 )p, 

(4.15) 

T ox -1 ox _pt ox -1 ox 
ot] 1 oro or) (mem) 52), 

_pt|%y-1 T| OX 7-1 ox 
Diy (M@M) = |p.+Dh Es mem) 5*),), 


As you can see, the discretization for two dimensions is more tedious than one 


dimension. 


The right-hand side of (3.6) is much simpler than the left-hand side. As a reminder, 
the right-hand side of (3.6) is 


2 (un) 
c fe" (Vp) ds. 


Since we are working with quadrilaterals, we integrate the four boundaries of 
each element, labeled by Qj. Figure 4.1 displays how the four boundaries of 
each element are numerically labeled. We pull the information from each side 


individually by using the following operators: 
L, =e) @1, L7 =1 ex, L3 =e @1, Ly =18e). (4.16) 


The four Ly operators (4.16) isolate the correct information from each element 
boundary. Additionally, the change of variables for the right-hand side produces a 
set of surface Jacobians found in either the € or 7 directions by 


i an 2 an 8 Ax 2 dy 2 : 
Sie=J (3") (3) = (= (3) j €=1,3, (4 17) 
2 2 2 2 
APR RPE ome 


By calculating the outward unit normal, n, for the sides of the element and incor- 
porating it in with the calculation of the numerical flux, the discretization of the 


right-hand side is 
7p" |LIMS iif 1 + LMS pf yp + L3MS jaf pg + LMS jaf pa], (4.19) 


where fy =1- (Vp). As a reminder from Section 4.1, $ ji Sj2, Sj3, and Sj4q are the 
matrix of surface Jacobians for each side of an element, as labeled in Figure 4.1. 
This same notation is used for all four sides of the flux term, fy listed in (4.19). 
The numerical flux for two dimensions is discussed in greater detail in Section 4.3. 
Finally, we combine the two sides of (3.6) to form the final discretization for the 


second variational equation to be 
dw 
J(M@M)—- + CT p =? (LIMS iif 9) + L3MS jf yp + L3MS af p3 + L{MS jaf p4)- (4-20) 


As a reminder, y does not appear in (4.20) because it must hold for all y, as 


discussed in Section 3.3.1. 


4.2.3 Third Two-Dimensional Discretization Equation 
The discretization process of (3.12) is similar to the discretization of (3.6) in the 


previous section. Starting with the left-hand side, 


op | 
VWV|—-@]dA, 
Lela 


we have eight terms from the product of the gradients 


voy (F- |- dé O (2 - \eoe dé O (F- oie 

ot Ox dé \ at Ox OE | Ax EN At ax an 
+33 (2 see eo (P- oo 
ax On\ a dx dE odxdn\oat ax an (4.21) 
dé O (Op déOW dé A [Op on ow ; 
+35 (Foor ale ol 
Pica 0 (op déOW dE O (Op on ow 
alee 0 eae sae |r 


oF 


Focusing on the first term of (4.21), we find the discrete form to be 


FY en 


-y"0t |S pment |p. Fo} (4.22) 


Expanding the discretization for all eight terms in (4.21) yields the full discrete 
form for the left-hand side of (3.12), 


dp dp 


T T sa 
“_qw)= ae & 4.2 
v'r(Z a) TS -p'Tw, (4.23) 
where T is equal to Equation (4.15). 
Let us focus on the right-hand side of (3.12), 
{. (Vw -n)(w* —w)ds, (4.24) 
aQ; 


and discretize side one of a single element (see Figure 4.1). Conducting a change 
of variables yields 


+1 
{ Vien foS dé, (4.25) 
-1 


where the fi, = (w* —w). For notational purposes, allow 


_ {dé on 

Dy = Ea + spy}, 
fog on 

Dy [SePe+ Su, | 


The discretization of (4.24) for side one becomes 


ip" |(LrDs)! mx +(L1Dy). ty] Sirf 
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where the flux, f 


similar calculation for the remaining sides gives 


wy iS from side one of an element as annotated in Figure 4.1. A 


ip" (DyL{ny+DyLi ny) MS fo, +" (Dy Lynx +DyLjny)MS pf .» 


(4.26) 
+" (DILIny + DIL ny) MS af 5+" (DELI n, + DELI my) MS jaf oa 


Equation (4.26) is the discrete form for the right side of (3.12). Combining (4.26) 
and (4.23) produces the final discrete approximation for the third variational equa- 
tion to be 

dp 


T—--—Tw= (D{Lin, a: DjL{n,)MS afer t (DIL; nyt DjL3n,)MS of on 


dt (4.27) 


Ty7T T7T T7T 
+ (DI LE ny + DIL ny) MS pf x3 + (DEL Ix + DELI ny) MS jaf oa. 


4.2.4 Combination of Two-Dimensional Discretization Equations 
As in one dimension, we have three discrete equations and our goal is to get the 
problem into the form of two equations and two vector unknowns. By combin- 
ing (4.10) and (4.27), we find 


d 
(T+1J(M@M) = =(T+1J(M@M) )\w+(D{Li ny +DyLi ny) MS iif 1 
+(D{ Lj nx + Dy Lj ny)MS pf .9 + (Dp Lynx +DyL3ny) MS af 3 (4.28) 
+(DILIn, + DIL n,)MS jaf 


wd 


where I is a matrix of ones since (4.10) produces a scalar. Additionally, for the sec- 
ond equation, (4.20) can be re-arranged to isolate deo for the final discrete equation 
of 


dw 
J(M@M)— =? ((L{MSji fp + L3MS jaf yp + L3MS jaf p3 + L4MS jaf p4) - Tp) (4.29) 


Equations (4.28) and (4.29) are the two-dimensional ordinary differential equations 
that we solve numerically using the RK54 scheme [11]. 
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4.3 Two-Dimensional Numerical Flux 

The numerical flux for two dimensions is conceptually similar to the numerical flux 
for one dimension, except in two dimensions we have to account for more degrees 
of freedom, boundaries, and the element normals. In Section 3.5, we accounted 
for the discontinuity that existed between the two different approximations at the 
boundary point. In this section, we account for the discontinuity that exists between 
two different approximations at multiple boundary points along the surface of two 
neighboring elements. Thus, the numerical flux in two dimensions is an extension 
of the numerical flux from Section 3.5. Figure 4.2 is an example depiction of a 


two-dimensional flux boundary for a quadrilateral element. 
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Figure 4.2: Example Two-Dimensional Flux Boundary on a Quadrilateral with N = 2. 


In two dimensions, the element to be considered is the (—) side and the neighboring 
element is the (+) side. This is the same for the (—) representing the element of 


focus’s gridpoints and the (+) representing the neighboring elements gridpoints. 
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In referencing Figure 4.2, the center element requires a flux to be computed with 
its neighbors at the boundary LGL points. Once again, Figure 4.2 is an example 
for discussion purposes; the actual elements are curvilinear quadrilaterals and the 
element normals at the LGL points are required in the numerical flux computation 
of (n-Vp)*. The formulation of the form of the central flux and upwind flux can be 


found in Section 3.5. The central flux equations for two dimensions are 
oe Repti 
w= 5 (@ +a"), (4.30) 
ods ss - + + 
n- (Vp) = 5 (n" -(Vp")-n* -(Vp")), (4.31) 


where (4.30) element normals are represented in (4.27). The element normals in two 
dimensions are outward of each element. For example, n~ is the normal directed 
away from the (—) boundary toward the neighboring (+) element and n* is directed 
toward the element of focus (i.e., the center element in Figure 4.2). Typically, the 
central flux is an average of both sides of an element boundary, but the reader may 
notice that (4.31) has a minus instead of a plus sign. The reason for this is that 
n* =—n", and thus the minus sign actually produces an average. As a reminder, 
the numerical flux is incorporated in the discretization equations by f, = -(Vp)’ 
from Section 4.2.2 and f,, = (w*—w) from Section 4.2.3. Similar to one dimension, 
the central flux is easy to understand and produces a stable algorithm, but it does 
not take into account the physical propagation of information like the upwind flux. 


The upwind flux equations for two dimensions are 


1 
a= 5 or +a*)—5(n*-(Vp*)+n7-(Vp")), (4.32) 
ine aes : 1 a 
n:(Vp) = 5 (n (Vp )—n*-(Vp")) +5 (@" -@ ). (4.33) 
Equations (4.32) and (4.33) are formulated by the decoupling of the wave equation 
into two one-way advection equations at the boundary points through the charac- 
teristic variables. This is the same concept from Section 3.5, but in two dimensions, 


there are three characteristic variables. However, the third characteristic variable 


is associated with a zero eigenvalue and does not propagate information across 
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element boundaries. As depicted in Figure 4.2, r; and rj are example characteristic 


variables propagating information across one elemental boundary. 


4.4 Two-Dimensional Discontinuous Galerkin Results 
As a test problem we consider a washer domain with curved elements as shown 
in Figure 4.3. Along the inner and outer radius of the washer, we impose the free 


surface boundary conditions of p = 0. We use the exact solution of 


p(x,y,t) = sin(nt— BO) Jp (r(x, y)), (4.34) 
Op 


f= w(x,y,t) =ncos(nt— BO) Jg(r(x,y)), (4.35) 


with the conversion equations 


r(x,y) = [x2 + y2, 


pe (ee 
O(x,y) = cos (5) 





for mapping between polar and Cartesian coordinates. Here, Jg is a Bessel function 
of the first kind with parameter f and n is an integer; we use 6 = 4 and n= 1. To 
ensure (4.34) is equal to zero at the inner and out radius for all time, we enforced the 
inner and outer radius of the washer grid to be the second and fourth roots of the 
Bessel function. For interested readers, the DG implementation for two dimensions 


can be found in Appendix C. 


We initially tested the implementation on a square grid of quadrilateral elements 
with periodic boundary conditions. After the successful implementation on a 
square grid, we built and tested the algorithm on a washer grid with curved 
elements using the exact solution of (4.34) and (4.35) evaluated at t = 0 as the initial 
condition. To enforce the free surface boundary condition of p = 0, we used 


wi =-w, (n*-(Vp*))=-(1 -(Vp_)) 


in the numerical flux routine for the points along the inner and outer radius of the 
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Figure 4.3: Example Two-Dimensional Washer Grid: Ne = 8 


washer. Once again, a visual depiction of an example washer grid with Ne = 8 
semi-curved elements is seen in Figure 4.3. Initially, testing the entire washer grid 
became too computationally expensive with the large number of elements. To 
counter this problem and still maintain integrity of the results, we evaluated one 
quarter of the washer by enforcing + periodic boundary conditions for the left and 
right boundaries through the bessel function. By letting 6 be any multiple of four, 
we enforced the + periodic conditions. For example, in referencing Figure 4.3, this 


means boundaries 1 and 3 from elements 1 and 2 are equal. 


Table 4.1 lists the different polynomial orders and elements tested using the initial 
condition with (4.28) and (4.29) and evaluated through a RK54 iterative method. 
Figure 4.4 shows the log of the error vs. log of the number of elements for N = 
2, 4, 6, 8 calculated using the global L? error norms. As expected for both w 
and p, Figure 4.4 displays increasing convergence rates as the polynomial order 
increases. Once again, the higher the order of the local approximation, the faster 


the convergence rates are due to the the global error being dependent on the 


45 


Table 4.1: Discontinuous Galerkin Tested Information 























Polynomial Orders (N) | Number of Elements (Ne) 
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Figure 4.4: Convergence Rates For N = 2, 4, 6,8 


polynomial order [1] (as discussed in Chapter 3 with (3.32)). 


Table 4.2 displays the convergence rates corresponding to Figure 4.4. As seen, 
the convergence rates increase with increasing polynomial order. Using the same 
space of functions for w and p, the two-dimensional results yielded convergence 
rates being near its associated polynomial order, with w being near N and p being 


N+ 5 or N+1. Once again, most DG methods are expected to be of the order N, 
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Table 4.2: Convergence Rates for N = 2, 4, 6, 8 















































Convergence Rates 
N/Ne_ | Ne=4tol6 | Ne=16 to 64 | Ne = 64 to 256 
N =2(a) 1.9240 2.0570 2.1274 
N=2(p) 1.9663 2.9243 2.8680 
N=4(a) 3.3438 4.2448 4.0822 
N=4(p) 4.2739 5.1808 4.4615 
N=6(a) 5.2521 6.0193 6.0944 
N = 6 (p) 6.6670 6.7889 6.4519 
N =8 (a) 6.9754 7.9190 8.0543 
N = 8 (p) 8.5157 8.7873 8.5069 

















N+ $, or N+1[1]. These results are different from the one-dimensional convergence 
rates, where w is N and p is N+2. Understanding this difference requires further 
analysis and is an area for future research. 


Table 4.3: Convergence Rate for N = 10, 12 
































Convergence Rates 
N/Ne Ne =4 to 16 | Ne = 16 to 64 
N = 10 (@) 8.7628 10.0287 
N = 10 (p) 10.2080 10.6994 
N = 12 (a) 10.6378 12.0772 
N = 12 (p) 12.0461 11.8315 











To use higher-order polynomials, we decreased the number of elements in order 
to avoid reaching machine precision with the first datapoint. Figure 4.5 shows 
the error and Table 4.3 gives the convergence rates for N= 10 and N=12. The 
convergence rates are still on the order of the polynomial, but the N = 12 case 
seems to be reaching machine precision. Both Figure 4.5 and Table 4.3 display how 
higher order polynomials can be used to approximate the solution with a much 


higher convergence rate. 


Using higher order polynomials, especially on complex geometries with curved 
elements, comes with a computational cost with respect to time for approximating 


the solution. Higher orders are more accurate, as seen in Figure 4.5, but the 
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Figure 4.5: Convergence Rates For N = 10, 12 


computational time increases since the matrices to be inverted are larger with 


higher polynomial orders and require a smaller time-step. In the code we use a 


time-step of At ~ 4. 
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CHAPTER 5: 
Energy Conservation 





Proving stability of a method can be accomplished through an energy analysis. 
The following energy analysis is similar to the method employed by Appel6 and 
Hagstrom [2]. By conducting the energy analysis on an element with neighboring 
elements, we are able to find out if energy is dissipated or conserved throughout 
the system. If energy is dissipated, then we know that all of the eigenvalues of 
the system have a negative real part and thus the ODE is stable. If energy (E) is 
conserved, then ap = 0 and all the eigenvalues are purely imaginary. If E; is the 


energy on an element k, then the total energy is 


Ne 
E= ye Ey, (5.1) 
k=1 
and we want to show 
Ne 

dE dE; 

— = ie A 

dt = dt ~ 0 Oe 


in order to prove that energy is dissipated throughout the system as time progresses. 


5.1 Basic Theory of Energy Conservation 


Let us look at the continuous energy equation on the domain x € [a,b] for one 


b 2 
= 1 5 1 (op 
E= | i” +3 (3 dx, (5.3) 


with p; = w, wt = CD xx, and c2 = = We assume c, A, and p are constants on the 


dimension: 





domain. As a reminder, for clarification purposes, ae = p; for a short notation form 
and this follows through for the remaining terms in Section 5.1. The first portion 


of (5.3) is the kinetic energy and the second portion is the potential energy for the 
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system. By taking the time derivative of (5.3) 


iE = [door dx (5.4) 
ioe Pe t pike y ; 


we can now begin to manipulate (5.4) to show (5.2). Specifically, we want to 
manipulate (5.4) in order to select our boundary conditions to ensure FE; < 0. Let us 
substitute w; = CDity and w = p; into (5.4) 


sie 1 
Ee= —CDxxPt + —Pxtpxdx. (5.5) 
a A p 
With c? = ao it follows that (5.5) becomes 


1 b 
care { PxxPt + prtpxdx, (5.6) 


and by conducting integration by parts on the last portion of (5.6) we find the 


4 b b 
Ey = | i Paxpidx + pxpt |? - { pups} (5.7) 


The two remaining integrals cancel and we have that 


following: 


1 1 
Ep = —pxp (P= [pro |p px lal, (5.8) 
t= OP ae al b —Px® a] 


where we can now choose the boundary conditions to ensure E; < 0. For example, 
with periodic boundary conditions, it should be obvious that (5.8) equals zero when 


a = b and energy is conserved. 


5.2 Two-Dimensional Energy Analysis 
Section 5.1 is a relatively simple example of the analysis of energy for the continu- 
ous one-dimensional system. Moving into two dimensions requires more tedious 


analysis and also depends on the conditions set within the system. For example, 
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when using an upwind flux, we want to show that energy is dissipated and when 
using a central flux, we want to show that energy is conserved. Let us initially 
investigate a general two-dimensional case before moving into the flux analysis. 
The energy equation for two dimensions is 


ae Pe ab Be Lerleoh 39 
B= { 50 55 (pi+p;)dA, (5.9) 
J 
where 
Op _ dé Op pipealidig 
Ox axoe ax On’ 
and 


ap) _ op zs “dp ap aeanap , Apdnaé dp | ap an ap 
par JE AEdxdAxOn OAndxdxdE An\A 


Ox dE Ax Ax AN 


2 2 
Moreover, (3) is formulated through the same method as (3) , except with 
respect to y. The change of variables of (5.9) yields 


wf, Le acters) a) 





dean | Ean CP ae, dE An | 9E On| Op (5.10) 
" Ox Ox * Oy dy On on | Ox dx * Oy dy dé 
3 é 
Sel) «C5 [eee 
and the discretization of (5.10) is 
ES =e T](MeM)a+ 5-9" Tp, (5.11) 
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where T is the same simplifying variable used in Equation (4.14) from Section 4.2.2. 
Taking the time derivative of (5.11) we have 


tis _ 1 (i 


T 
1 of dw 


which can be simplified to 


és 1 rragoa(4@\alora(® 


Upon inspection of (5.12), we can substitute (4.27) and (4.29) in for J(M@M) (42) 


and T(£) from Chapter 4 to produce 


dEk Creer T F i 
Ghee [(LTMS 1 fp + LMS 2 f pp + L3MS jaf y3 + L{MS ja ian -Tp| 
1 
+ sas [Tw +(D{Liny +DyLi ny) MS if.) + (Dy Lynx +DyL3 ny) MS pf 9-19) 
T7T T7T T7T T7T 
+ (Di L3nx + Dy L3ny)MSjsf 43+ (Dy Linx + Dy Lj ny) MS jaf 4]: 
For notational purposes, let 
_(7T ‘g T Tape: 
F, = (LTMS jf, + LZMS of 9p + L3MS af 3 + LIMS sf ya), 
F,, =(D{L{ ny +DyLi ny) MS iif.) + (Di Lynx +Dy Lj ny) MS pf 9 
+(D{L3ny+ DjL3ny)MS jafog+(DrLinx+ DjLiny)MS a ne 


Therefore, Equation (5.13), with the substitution of c = 3 becomes 


a 5° (F, 7 Tp) + pl (re +F,), (5.14) 
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and since w! Tp = p' Tw, (5.14) is further simplified to 


7 5 Fy +p! Fa). (5.15) 


As you can see, (5.15) is dependent on the flux and boundary information from 
each element, as included in F, and F,,. Knowing that 
Ww, = w'Li, 
(ny . (vp,)). = p! (D{Lin, + Din) , 
for €=1,2,3,4, (5.15) becomes 


dE, 1; ¢ i £ T 
an = 5 erMSaf +0 MS j2f 9 + @3MSi3f3 + @4MS jaf 


+(m- (vp,)) M Sitfor + (m2: (vp,)) MSix ie (5.16) 
T iy 
” (13 -(Vps)) MS jsf 3+ (1a: (Vp,)) MS af os|: 


Recall that f, =n-(Vp)' and fi, = (w* — w) from Section 4.2.2 and Section 4.2.3, thus 


we can consolidate the boundary information for each side of an element 


Fs = 1 cof (n-(Ve)) +(0-(Vp,))" MS (o — er) 


+ (oIMS (12 : (Vp,)) a (n2 ; ( 
N00 


+ (ofMS j (114 . (Vp,)) + (114 (Vp, ) MS j4(@* — w))| : 


T 
Vp») MS j2 (w" - @2)] 5.17) 
+ (ofS j3(1-(Vp, (Vps)): MS} (co" - ws) 
) T 


Equation (5.17) is in the final form that we use to conduct the stability analysis. 
Summing (5.17) for all elements, we can analyze the energy for the whole domain. 
In what follows, we prove stability with the boundary conditions imposed from 


Chapter 4 for the central flux and the upwind flux. 
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5.2.1 Central Flux Analysis 

In order to simplify the analysis, we focus on a single side of an element with its 
adjacent neighboring element. This concept is depicted in Figure 3.2 in Chapter 3. 
For notation purposes, the (—) side is the element of focus and the (+) side is the 
neighboring element. Additionally, w~ represents the solution at the degrees of 
freedom on the associated side of the quadrilateral element as depicted in Figure 4.2. 
This same notation applies for (5.19) and (5.20). The equations for the central flux 


for a generic side without boundary conditions are 


w= s(o +0"), (5.18) 
1 

a -(Vp')= 5 [0 -(Vp)—n"-(Vp")], (5.19) 
1 

n*-(Vp") = 5[0' -(Vp")—n (Vp). (5.20) 


Since we are working with one side of a quadrilateral element, we add together 
the terms from one edge of the element with the terms from the neighbor’s corre- 





sponding edge 
TE = [(*) MSj(n*-(Wpy')+(n" (9p) MS;(w" - 0°) 
de NY pee LE NE ier Xen (5.21) 
+|(w~)' MS;(n~- (Vp)') + (17 -(Vp-))' MS;(w* - @)]. 
Substituting in (5.18), (5.19), and (5.20) for the flux portions in (5.21) yields 
dE 1 +yT + + + +)\yT + + +)\\T - 
—E <5 [(@*)MS\(n* (Vp) — (nt -(Vp*)) MS jco* + (n* (Vp) MS jo 
_ (5.22) 


—(w-)' MS;(n* -(Vp*)) +(@)' MS; (n--(Vp-))- (0 - (Vp) MS jw 
+(n~-(Vp~))' MS;w* -(w*) MS;(n" - (Vp-))]. 


d = 
Inspecting (5.22), all of the terms cancel out to show a = 0 for one element 
boundary. This can be expanded to all four boundaries of the quadrilateral. Thus, 


energy is conserved and the ODE is stable when using the central flux. 
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5.2.2. Upwind Flux Analysis 

As discussed at the end of Section 3.5, the upwind flux consists of the central flux 
plus an upwinding portion. This is evident in (3.28) and (3.29). Therefore, since we 
know that energy is conserved (i.e., ae = 0) for the central flux portion, we only 
investigate the upwinding portion. The equations for the upwinding section are 


wi, =—S[n*-(Vp +0 (Vp), (5.23) 
(n* (Vp), = 5-(@" -"), (6.28) 
(1 -(Vp)'), = 5-(@* 0°), (5.25) 


which is substituted into the elemental boundary equation (5.26) for the upwinding 


portion, 
dE* 
al =|(@*)' MS; (n* -(Vp)’), + (a* -(Vp*))" MS 0%, (5.26) 
+|(@~)' MS;(n" -(Vp)"), + 0 «(Vp ))' MS a]. 
The substitution of (5.23), (5.24), and (5.25) into (5.26) yields 
dE* 
(=! =5- [(w*)' MS;(w" - w*) + (w")' MS;(w* -w)| 
: (5.27) 


c a 2 = 
— 5 [(0*-(Vp*))' MS; [n* -(Vp*) +0 - (Vp) 
= eT 7 _ 
+ (0 -(Vp"))' MS;[n* -(Vp*)+n--(Vp)]]. 
Factoring out a negative from the first portion of (5.27) makes (w7 — w*) become 


(w* —w). Therefore, (5.27) can be re-written to 


dE* 1 es ay + = 
(Gt) =- sel -w") MS;(w* -w)| (5.28) 


— 5 |" (pt) +a (Wp) MS;(a* (Vp) + - (Vp) ]- 
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Inspection of (5.28) displays two terms: 


[(w* - w-)' MS; (w* — w)| ’ 
[(n* -(Vp*) +n" -(Vp-))' MS; (n* -(Vp*) +0 -(Vp-))| 


that are both positive definite. MS; in both sections is always positive; therefore, re- 
gardless if (w* —w7) or (n*-(Vp*)+n7 -(Vp7)) produce negative or positive values, 
the square of each is always be positive. Furthermore, with the negative in front 
of both sections of (5.28), (=) <0. The combination of (=) with aa from the 
central flux section (Section 5.2.1) yields a negative value per face. The summation 
of all the elements in the domain shows Equation (5.2) to be true. Therefore, this 
method is stable for both the central flux and upwind flux on complex geometries 


with curved elements. 


5.2.3 Boundary Condition Analysis 

In Chapter 4, we tested the method on a washer mesh with curved elements and 
the boundary conditions of p = 0 on the inner and outer radius of the washer. In 
order to implement the boundary conditions, the following constraints were set at 


the inner and outer boundaries of the washer 
w =-w, (n'-(Vp"))=—-(n -(Vp’)). (5.29) 


Additionally, we are only concerned with the (—) portion of (5.21) for the central 
flux and the (—) portion of (5.26) for the upwind flux since there are no neighboring 


elements ((+) portions) at the inner and outer radius of the washer. 


Boundary Conditions with Central Flux 
Once again, substituting in the central flux Equations (5.18) and (5.19) for the (—) 
portion of (5.21) yields 


dE- 
SE = (7) MS)(5 [> (Vp) —n*-(Vp*)]) + O-(Wp))" MS;(5(w*-@)), 
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which becomes 


GE: 
FE = oJ! MS\(£ 1200 (Wp -))]) +r (wp MS;(4(-207)) 6.30) 


after the substitution of the boundary conditions listed in (5.29). Simplifying (5.30) 


produces 





ee T 2 
SE = (wo) MSj(n"-(Vp")) - (8 “(Wp ))MSj(w") =0, (6.31) 
and energy is conserved with the boundary conditions of p = 0 using a central flux. 


Boundary Conditions with Upwind Flux 

For the upwind flux boundary condition analysis, we use the same process as 
above, except we use (5.23) and (5.25) for substitution into the (—) portion of (5.26). 
This yields 





(FE) =o" Ms(L(w*-« )) +0" -(vpyMs,(—S[0" (wp) +9" (wp) 


which becomes 








(SE) = oryhas)( $20 )) + (--(vp-)™MS;(-S[-n-- (Vp) +n7-(Vp)]) 
and is simplified to 


(i 


| = -* [(w")' MS;w" | <0. (5.32) 


Once again, the combination of (5.31) and (5.32) produces a negative result for the 


upwind flux boundary condition analysis. 


The global energy dissipation rate is the sum of the energy dissipation rates at all of 
the faces. Therefore, for both the central and upwind fluxes, the method discussed 


in Chapter 4 is stable. 


D/- 
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CHAPTER 6: 
Conclusion 





The complexities of solving PDEs through Finite Difference Methods, Finite Volume 
Methods, and Finite Element Methods is and will continue to be an area of ongoing 
research. In this paper, we explored a new Discontinuous Galerkin method for 
approximating a second order wave equation with curved elements in complex 
geometries. The beginning two chapters discuss the tools needed for constructing 
the DG method, while Chapter 3 and Chapter 4 tested the method in one and two 
dimensions. Chapter 5 proved energy stability of the method through an energy 


analysis. 


In Chapter 3, convergence rates on a one-dimensional grid for both low and high 
polynomial orders yielded rates near the associated polynomial order. The results 
yielded convergence rates with w being on the order of N and p being on the order 
of N +2, which is higher than most DG methods that yield results on the order of 
N,N+ 5, and N+1. Additionally, using higher polynomial orders required fewer 


elements for a given error level. 


In two dimensions, we tested the problem on a washer grid with curved elements. 
In this case, the inner and outer radius had the free surface boundary condition 
of p = 0 implemented through the numerical flux routine. The two-dimensional 
results yielded convergence rates different than the one-dimensional method with 
w and p being near its associated polynomial order N or N+ 5. Similar to the one- 
dimensional results, testing high polynomial orders achieved machine precision at 


a faster rate with less elements. 


Chapter 5 explored the energy dissipation of the method, with and without the 
imposed boundary conditions of p = 0, for the central flux and upwind flux. In all 
cases, the method either conserved energy (2 = 0) or dissipated energy (4 < 0). 
Therefore, the eigenvalues of the system have a negative real part proving that the 
ODE is stable and the method is viable. 


oo 


6.1 


Future Work 


Future research into this Discontinuous Galerkin method can be conducted on a 


variety of topics. These topics include, but are not limited to, the following: 


Accuracy: As discussed in Chapter 3 and Chapter 4, the convergence rates 
listed are observations; proving these results could lead to further insight and 
improvement. 

Efficiency: Is it possible to template the curvilinear elements so that the mass 
matrix is the same for all elements by modifying the approximation space for 
each element? The benefit of doing this is the ability to store one mass matrix 
for all elements [12, 13]. 

Coupling: Development of the numerical coupling procedures for elastic— 
acoustic interfaces within this second order form [14]. 

Form: Relationship between this second order formulation and the standard 
first order formulation of the acoustic wave equation. 

Extension: Consider this method for other equations, such as the Einstein 
equations governing black hole dynamics where there are 10 equations in the 


second order form versus 50 equations in the first order form [15]. 
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APPENDIX A: 
Interpolation and Integration 





The following are the main MATLAB codes for Interpolation and Integration from 
Chapter 2. 


A.1 Interpolation 








%This is the main driver for Interpolation. 


3/%Written by Benjamin Davis 


Yo Department of Applied Mathematics 
5|%o Naval Postgraduate School 
% Monterey , CA 93943-5216 
%o 
%*Synopsis: Conducting Interpolation using LGL points. 
5 a a a a % 


%lnterpolate a known function f(x) using legendre-gauss-labatto points. 
clear 
%Input Nth order interpolation (N) and number of evaluation points (k) 


3IN = 40; 


k = 50; 


5)%Initialization 


errLInum = zeros(N,1); 


7}errLlden = zeros(N,1); 


errL1 = zeros(N,1); 


ojerrL2num = zeros(N,1) ; 


errL2den = zeros(N,1); 
errL2 = zeros(N,1); 


errinfabsn = zeros(k,1); 
s)errinfabsd = zeros(k,1); 

errinf = zeros(N,1); 

for n = 1:N 


*Setting up Grids 
x = legendre_gauss_lobatto(n+1); 
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ou 





z = linspace(-1,1,50); 

%Set up data for Lagrange Li(Xk) 

f = exp(-4#x.42); 

%Constuct Lagrange 

[L,dL] = lagrange_basis(x,z); 

%Evaluation 

Pn = f « L; 

%Error Analysis 

fex = exp(-4+z.%2); 

for i = 1:k 
errLlnum(n) = abs(Pn(i)-fex(i)) + errLinum(n) ; 
errLlden(n) = abs(fex(i)) + errLlden(n) ; 
errL2num(n) = (Pn(i)-fex(i))%2 + errL2num(n) ; 
errL2den(n) = (fex(i))42 + errL2den(n) ; 
errinfabsn(i) = abs(Pn(i)-fex(i)); 

errinfabsd(i) = abs(fex(i)); 

end 


errL1(n) errLinum(n)/errL1iden(n) ; 
errL2(n) = sqrt(errL2num(n)/errL2den(n) ) ; 
errinf(n) = max(errinfabsn) /max(errinfabsd) ; 


SYS/S/S/SISYS/S15/5/8/81513/5/5/5/8/8/81515/5/8/8/81515/5/5/5/8/818185/5/5/8/8/81515/5/5/5/81515/5/5/5/0 


1 2 3 4 5 6 7 8 
[1 0 0;0 0 0;0 0 0; 010; 000; 00 0; 00 0; 10 O;... 
00 0;0 0 0;0 0 0; 00 0; 00 0; 00 0; 00 0; 0.5 0.3 0.2]; 
9 10 1 12 13 14 15 16 


0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0, 
0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0 


Yplot Pn vs. Actual Function 

hold on 

plot(x,f,'«' ) 

title ('Lagendre Gauss Lobatto Approximation ') 
xlabel('x') 

ylabel('f(x)') 

axis([-1 1 -0.1 1]) 

legend('N = 2', 'N= 4','N= 8', 'N= 16', ‘Exact') 


%Plot Error Norms 


figure 
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N = 1:40; 

69 semilogy(N,errL1 ,N,errL2 ,N, errinf ) 

title ('Legendre Gauss Lobatto Interpolation Error') 
71 xlabel('N') 

ylabel('Error Norm') 

73 legend('L1 Norm','L2 Norm', 'Inf Norm') 











A.2 Integration 





%This is the main driver for Integration. 


w 


%Written by Benjamin Davis 


Yo Department of Applied Mathematics 
5|% Naval Postgraduate School 
% Monterey, CA 93943-5216 
7|% 
*Synopsis: Conducting Integration using LGL points. 
2 % ee gg tagP me ist ey SOM a te rte Ae MS VEST an St SOOT ot ee bane EPS IE (ES See epee Aten Oe MSP (Ean) ange a Emule ie ew ran at eT EPC Part Sar NEE it ge art elgie T % 


%lntegration: Using the interpolation functions, sampling the 


B 


%basis functions at LGL and LG integration points to perform the Gauss 
Yquadrature of the given equation from the project set. 
clear 


wo 


% Initialization 


a 





N = 19; 
7 
for n = 1:N 
9 %Quadratrue points and weights 
[x,w] = legendre_gauss_lobatto(n+1); 
21 %Given Equation for evaluation 
f = exp(-4#x.42); 
23 %Constuct Lagrange matrix and differentiation matrix 
[L,dL] = lagrange_basis(x,x); 
25 %Evaluation 
Pn = (f+#L)+ w'; 
27 %Evaluation of Error 
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31| end 





exact = (sqrt(pi)/2)serf(2); 
errL1(n) = abs(Pn - exact)/abs(exact) ; 
errL2(n) = sqrt((Pn - exact) %2/(exact)%2); 


%Plot Error 
N = 1:19; 
semilogy(N, errL2) 


title ('Legendre Gauss Lobatto Integration Error’) 
xlabel('N') 


ylabel('Error') 
legend('L2 Norm') 
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APPENDIX B: 
One-Dimensional Discontinuous Galerkin 





The following are the main MATLAB codes for the one-dimensional problem from 
Chapter 3. 








%This is the Driver function using a Discontinuous Galerkin method for 
approximating a 2nd Order Acoustic Wave equation in one dimension 
Ywith an Upwind Flux. 


% 

%Written by Benjamin Davis Created: October 2014 
%o Department of Applied Mathematics 

% Naval Postgraduate School 

% Monterey, CA 93943-5216 

% 


*%Synopsis: Discontinuous Galerkin Method for wave equation in second 
*order form using inexact integration an Upwind flux. The outputs 
Yare currently four plots: Convergence rates for w and P and plots 
%of the numerical and exact solutions. 


clear 

%Initial Inputs 

N= 4; %Polynomial Order 
n = 3; 

dtscale = 1/4; 

c= 13 

Z =A 

t_final = 0.58; 

ngl = N+1; 


for Ne = 2.4(1:4) 
fprintf('Number of Elements %4d with polynomial order %2d\n' ,Ne,N) 
Np = Nesx(ngl); 
%Interpolation and Integration Points 
[psx,w] = legendre_gauss_lobatto(ngl) ; 
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40 


44 


46 


48 


Oo 
np 
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ea 





%Constuct Lagrange matrix and differentiation matrix 
[L,dL] = lagrange_basis(psx,psx) ; 

%Construct Mass and Differentiation Matrices 

[M,D] = mass_diff1D(L,dL,w) ; 

%Construct Global Mass and Differentiation Matrices 


[coord,intma] = create_grid(ngl ,Ne, psx) ; 
dx = coord (N+1)-coord (1) ; 

M = M.«dx/2; 

D = M\D; 

Dh =D"; 


grad2 = DhsM:sD; 
one = ones(ngl,ngl) ; 
onesM = one:M; 


M1 = grad2 + onesM; 


e0 = sparse(1,1,1,ngl,1); 


en sparse(ngl,1,1,ngl,1); 
%*Building the Initial Condition 
for e = 1:Ne 
for i = 1:N+1 
I = intma(e,i); 
x = coord(I); 
qO(e,i,1) = nspixsin(nspi+x); %dp/dt=w at t=0 
qO(e,i,2) = 0; %P at t=0 
end 
end 
ql = qO(:,:,:); 
%lime Step Calculation 
dt = dtscale+(1/Ne) /(N‘%2); 
Nt = round(t_final/dt); 
dt t_final/Nt; 
%Periodic Boundary Conditions 
[sidep] = sideper (Ne); 
%Facemap 


[fmp] = facemap (Ne) ; 


0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/9, 
0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0 
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N 


80 


84 


88 


90 


00 


08 





PIS/S/S/SSIS/S/S/5/51 5151515151315 /8/5/5/51515181515/5/5/51551518/5/3/5/5/5/5/5155181518/5/5/5155181815/3/5/5/5/5/5(5/818/5/8/5/5/0 
YCompute RK Time-Integration Coefficients 
kstages = 4; %=1 is RK1; =2 is RK2; = 3 is RK3; and = 4 is RK4 
if(kstages < 4) 
[a0,al,beta] = compute_ti_coefficients (kstages) ; 


for k = 1:Nt 
for a = 1:kstages 
[wf,dpf] = upwindflux(Ne,q1,D,ngl,sidep) ; 
[R] = RHSDG(en,e0,Dh,D,M,q1,Ne,MI1,wf,fmp,dpf,c) ; 


qm=a0(a)+q0 + al(a)+ql + dtsbeta(a) +R; 


ql = qm; 
end 
q0 = qm; 
end 
elseif(kstages == 4) 
a = [0,1/2,1/2,1]; 
b = [1/6,1/3,1/3 ,1/6]; 
R= 0; 
for k = 1:Nt 
qm = q0; 
for s = 1:4 
qs = q0+dtsa(s)+R; 
[wf,dpf] = upwindflux(Ne,qs,D,ngl,sidep) ; 
[R] = RHSDG(en,e0,Dh,D,M, qs ,Ne,M1,wf,fmp,dpf,c) ; 
qm = qm + dt«b(s)«R; 
end 
q0 = qm, 
end 
end 


®PYHSSSLSSS SSS SSS SSSSS SSS SSS SSS SSS SSS SS SS SSS SSS SV SSSS SSS SS SS SSS SSIS 
0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0 


®PBVASYS SYS SASS SSS SSSS SSS SSSS SSS SSS SASS SS SS SSS SSSSSSSSSSSSSS SSS SSSSSSS 
0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0 
é ; 
%Plotting Purposes 
for e = 1:Ne 


for i = 1:N+1 
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120 


122 


126 


128 


144 





end 


I = intma(e,i); 
x = coord(I); 
qex(e,i,1) = nspixsin(nspi+x)+cos(nspi+t_final); 
qex(e,i,2) = sin(n+pi+x)*sin(nspist_final); 
end 

end 

%Plotting purposes 

m=1; 

for e = 1:Ne 

for i = 1:N+1 

w(m) = qm(e,i,1); 
wexact(m) = qex(e,i,1); 


p(m) = qm(e,i,2); 
pexact(m) = qex(e,i,2); 


I = intma(e,i); 
xp(m) = coord(I); 
mem+1; 
end 
end 
%Building the Error Plot Information 
L2errw = 0; 
L2errp = 0; 
for elm = 1:Ne 
pnts = (elm-1)+(N+1) + (1:(N+1)); 
dw = w(pnts)-wexact(pnts) ; 
dp = p(pnts)-pexact(pnts) ; 
L2errw = dw«Medw'+ L2errw ; 
L2errp = dpsMedp'+L2errp; 


end 

L2err(1,z) = Np; 
L2err(2,z) = sqrt(L2errw) ; 
L2err(3,z) = sqrt(L2errp) ; 
L2err(4,z) = Ne; 

z= z+; 
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8| % YASYSISYS/S/3/5/8181515/5/5/8/5/81515/5/5/8/81518185/5/5/8/5181515/5/8/5/8151513/5/5/8/5/81515/5/5/8/8/51515/5/5/5/8/5/o 


% 1 2 3 4 5 6 7 8 
c = [1 0 0;0 1 0;0 0 0; 00 0; 00 0; 00 1; 0 0 0; 10 0O;... 
00 0;0 0 0;0 0 0; 00 0; 00 0; 00 0; 00 0; 0.5 0.3 0.2]; 


52| % 9 10 ate 12 13 14 15 16 


Yo. ASISLSLS/S/S/SISISIS/S/S/S/S/S/SIS/S/S/S/S/SISISLS/S/S/S/S/S/S/S/S/S/S/5/81 S181 5/5/5/5/5/S1 S/S S/S/S/S/S/S1S/S/S/S/S/S/SISISLo 
PIS/SISISISISISISISISISIN Convergence P10 t/S/S/S/S/SISISISISISISIS/SIS/S/S/5/5/S/5/S/S/S131S1o 
figure (1) 
subplot(1,2,1); %hold on 
loglog(L2err (4,:) ,L2err(2,:),'*-','color',c(N,:)) 
title (‘Convergence Plot for w') 
xlabel('Ne - Number of Elements ') 
ylabel('L2 w') 
axis tight 
legend('N = 4', 'N = 6','N = 8','N = 16') 
PIISIISISISISIIIIIE Convergence P10 t/s/S/S/S/S{S/SISISISISISISIS/S/S/5/5/S/5/S/S/S/81S1o 
subplot(1,2,2); %hold on 
loglog(L2err (4,:) ,L2err(3,:),'*-','color',c(N,:)) 
title (‘Convergence Plot for p') 
xlabel('Ne - Number of Elements ') 
ylabel('L2 p') 
axis tight 
legend('N = 4', 'N= 6', 'N= 8','N= 16') 
% WISSISISISISISISS/S/S/S/0P lot Numerical Solution P(x, t )YWYSS/S/SISISISISISo 
figure (2) 
plot(xp,p,xp,pexact,'»') 
title ('Numerical Solution for p(x,t)') 
xlabel('x') 
ylabel('t') 
legend( 'Numerical','Exact') 
%o VYSISNSTSISINSVSISIS/S/0R lot Numerical Solution w(x, t )¥S/S/S/S/SISISISIS/S/o 
figure (3) 
plot (xp,w,xp,wexact, ‘x ') 
title (‘Numerical Solution for w(x,t)') 
xlabel('x') 
ylabel('t') 
legend( 'Numerical','Exact') 


6jrate_w = (log(L2err(2,1:end-1))-log(L2err(2,2:end))) ./ 
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r= 





(log (L2err(1,2:end))-log(L2err(1,1:end-1))) 


s|rate_p = (log(L2err(3,1:end-1))-log(L2err(3,2:end))) ./ 


(log (L2err(1,2:end))-log(L2err(1,1:end-1))); 


fprintf('Convergence rate for w: ') 


2| disp (rate_w) 


fprintf( ‘Convergence rate for p: ') 
disp (rate_p) 








%Function for the Right Hand Side of Discretized equations for 


3}%one dimension DG from Chapter 3. 


%Written by Benjamin Davis Created: November 2014 


5|% Department of Applied Mathematics 

Yo Naval Postgraduate School 

71% Monterey , CA 93943-5216 

% rhe eles on geet ee et ae ema ey eee tw dee Menem, ae ees meer ee ome oe See ear Leretael a! aha ke Pam er Se me a ee en ee tar eye ee er et ae on! a Gar ey ee wages oe % 





function [R] = RHSDG(en,e0,Dh,D,M,qn,Ne,M1,wf,fmp,dpf,c) 
for e = 1:Ne 
R(e,: ,2) Mixqn(e,:,1) '; 
R(e,: ,2) R(e,:,2)' + Dhs(enswf(fmp(2,e))-e0swf(fmp(1,e))); 
R(e,:,2) = R(e,:,2)' - Dhs(ensen'sqn(e,:,1) '-e0+e0'sqn(e,:,1) ') 


R(e,:,1) = -c42sDh+sM+Dsqn(e,:,2) '; 
R(e,:,1) R(e,:,1)' + c42s(en+dpf(fmp(2,e))-e0sdpf(fmp(1,e))); 


A 
any 
Il 


MI\R(:,:,2) '; 
R1 = R1'; 
R(:,:,2) = R1; 

R2 =M\R(:,:,1) '; 


23} end 
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%Function for the Center Flux for the one-dimensional DG problem from 
%Chapter 3 
%Written by Benjamin Davis Created: October 2014 


5|%o Department of Applied Mathematics 

% Naval Postgraduate School 

7|% Monterey , CA 93943-5216 

% tear A art es etme eee vene Mae a Pee epee aes ee ed epee Ne Pan Aa Wee et ras ora rae eh ce oes ee melee nme” Seer Ane ee an rae eS nae een eg ee” Se ney en Pe Ege ee rs oe aeey eee age ae epee ee ear ta % 





w 





function [wf,dpf] = centerflux(Ne,qn,D,ngl,sidep) 
for s = 1:Ne 
Ls = sidep(1,s); 
Rs = sidep(2,s); 


qLw = qn(Ls,ngl,1); 
qRw = qn(Rs,1,1); 
dpL = Dsqn(Ls,:,2) '; 
dpL = dpL(ngl); 

dpR = Dsqn(Rs,: ,2) '; 
dpR = dpR(1); 


wf(s,:) = (1/2) +(qLw + qRw) ; 
dpf(s,:) = (1/2)*(dpL + dpR); 
end 


s}end 








%Upwind Flux function for one dimension DG. Used in one dimension driver 
%for 2nd Order Acoustic wave equation. 
%Written by Benjamin Davis Created: October 2014 


5] %o Department of Applied Mathematics 

Yo Naval Postgraduate School 

71% Monterey, CA 93943-5216 

% OR ee ie aay Oe ee he gy ee ee ee ee Ee Gey ee ee ee ee ee cee eee ee ene % 





function [wf,dpf] = upwindflux(Ne,qn,D,ngl,sidep) 
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nN 


» 


a 








c= 1; 
for s = 1:Ne 
Ls = sidep(1,s); 
Rs = sidep(2,s); 


qLw = qn(Ls,ngl,1); 
qRw = qn(Rs,1,1); 


dpL = Dsqn(Ls,: ,2) '; 
dpL = dpL(ng]l); 

dpR = D+eqn(Rs,:,2) '; 
dpR = dpR(1); 


wf(s,:) = (1/2) «(qLw + qRw) + c/2+(dpR - dpL); 
dpf(s,:) = (1/2) *«(dpL + dpR) + 1/(2+c)*(qRw - qLw); 
end 
end 











yCode given to Professor's F.X. Giraldo MA4245 Class July 2014 
%Used by Ben Davis 

%This code computes the Legendre-Gauss-Lobatto points and weights 
Ywhich are the roots of the Lobatto Polynomials. 


(|\%Written by F.X. Giraldo on 4/2000 


%o Department of Applied Mathematics 

% Naval Postgraduate School 

% Monterey , CA 93943-5216 
RR Yo 
function [xgl,wgl] = legendre_gauss_lobatto(P) 

p=P-1; 


ph=floor( (p+1)/2 ); 


for i=1:ph 
x=cos( (2+i-1)#pi/(2sp+1) ); 
for k=1:20 
[LO,LO_1,L0_2]=legendre_poly (p,x) ; 


72 








6 


34 


w 
oe 





%Get new Newton Iteration 
dx =-(1-x%2)sL0_1/(-2+x«LO_1 + (1-x‘%2)sL0_2); 
X=x+dx ; 
if (abs(dx) < 1.0e-20) 
break 

end 

end 

xgl(p+2-i)=x; 

wegl(p+2-i) =2/(p*(p+1)+«L0%2); 

end 


30/%Check for Zero Root 


if (p+1 ~= 2«ph) 
x=0; 
[LO,LO_1,LO_2]=legendre_poly(p,x); 
xgl(ph+1)=x; 
wegl(ph+1) =2/(p*(p+1)+L0%2) ; 
end 
%Find remainder of roots via symmetry 
for i=1:ph 
xgl(i)=-xgl(p+2-i); 
wel(i)=+twel(pt+2-i); 
end 








%Function for building the Lagrange Polynomials. 


3|%Written by Benjamin Davis in MA4245 Created: July 2014 in MA4245 


% Department of Applied Mathematics 
5|% Naval Postgraduate School 
% Monterey, CA 93943-5216 
va % eS a2 Se a aise Se SSS SS eee a eS See ee ee a a a a a a ee aS ae Sa a ea ee a % 


oO 


B 





function [L,dL] = lagrange_basis(x,z) 

%Nth order interpolation 

n = length(x); 

%Length of the equally spaced grid for k = 1:50 
h = length(z); 
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a 





%Initialize the Lagrange Matrix 
L = ones(n,h); 

dL = zeros(n,h); 

%Computation for Lagrange Matrix 


for k = 1:h 
for i= 1I:n 
for j = 1:n 
dl = 1; 


if j ~= i% If j does not equal i 
%Equation for the Lagrange Polynomial 
L(i,k) = (z(k)-x(j)) -/(xGi)-x(j)) + LG,_k); 
for 1 = 1:n 
if (1 ~= i) & (1 ~= j) 
dl = di«(z(k)-x(1)) ./(x(i)-x(1)); 











end 
end 
dL(i,k) = dL(i,k) + dl/(x(i)-x(j)); 
end 
end 
end 
end 
end 
a Yo 
%Code given by Professor F.X. Giraldo to MA4245 class. 
%Used by Ben Davis 
%This code computes the Legendre Polynomials and its 1st and 2nd 
%derivatives 
%Written by F.X. Giraldo on 4/2000 
% Department of Applied Mathematics 
%o Naval Postgraduate School 
% Monterey , CA 93943-5216 
% 
%This code was written by Professor Giraldo and given to his MA4245 
%Galerkin Methods class in July 2014. 
Yor a2 ccc cc esc c cee cern cess cree sce e nc se scree ese rc eee cece s secs nse ssccces Yo 








function [LO,L0_1,L0_2] = legendre_poly(p,x) 


S 


5} L1=0;L1_1=0;L1_2=0; 
LO=1;L0_1=0;L0_2=0; 


for i=l:p 

20 L2=L1;L2_1=L1_1;L2_2=L1_2; 
L1=L0;L1_1=L0_1;L1_2=L0_2; 

2 a=(2si-1)/i; 

b=(i-1)/i; 

24 LO=axx*Ll1 - b«L2; 

LO_l=a»(L1 + x+#L1_1) - beL2_1; 

26 LO_2=a%(2+L1_1 + x+L1_2) - bxL2_2; 
end 




















%Function for building the Mass and Differentiation matrices for 
3|%One-Dimension. Used in Thesis 1D Upwind code. 
%Written by Benjamin Davis Created: July 2014 


5|%o Department of Applied Mathematics 

Yo Naval Postgraduate School 

71% Monterey , CA 93943-5216 

% ee ee gee ge ye ee gee a ae ee mm ee ee ee ee ee ee ee pee me eee Oe Yet % 


function [M,D] = mass_diff1D(L,dL,w) 
n = length(L(:,1)); 
iM = zeros(n,n); 


D = zeros(n,n); 


for i = 1:n 

5 for j = 1:n 
M(i ,j) 
7 D(i,j) 
end 


(LG ,:) «LG ,:) )ew(1,:)'); 
((L(i,:) .*dL(j ,:) )sw(1,:)'); 





9 end 
end 
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yCode given to Professor's F.X. Giraldo's MA4245 class 
%Used by Ben Davis 

4/%This function computes the LGL grid and elements. 
%Written by F.X. Giraldo on 10/2003 


6| Yo Department of Applied Mathematics 

Yo Naval Postgraduate School 
8} % Monterey, CA 93943-5216 

%G- -- ----- ~~ 2-2 2 oe 2 2 --- % 
o} function [coord,intma] = create_grid(ngl,nelem, xgl) 


%Set some constants 
2}xmin=-1; 

xmax=+1; 

4| dx=(xmax-xmin) /nelem ; 
%Generate Grid Points 
6 ip =1; 

coord (1)=xmin; 





s|}for i=1:nelem 
x0=xmin + (i-1)+dx; 














20 intma(i,1)=ip; 
for j =2:ngl 
2 ip=ip + 1; 
coord(ip)=( xgl(j)+1 )*dx/2 + x0; 
24 intma(i,j)=ip; 
end 
2} end 
%G-------- ~~~ 22-2 2 2 2 eee % 


2}% Function for implementing periodic boundary conditions for Thesis 1D 
% Upwind code. 
4J%Written by Benjamin Davis Created: November 2014 
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%o Department of Applied Mathematics 











%o Naval Postgraduate School 
% Monterey , CA 93943-5216 
Re Yo 
function [sidep] = sideper (Ne) 
sidep = zeros(2,Ne); 
for e = 1:Ne 
sidep(l,e) = e-1; 
sidep(2,e) =e; 
if e==1 
sidep(1,e) = Ne; 
end 
end 
end 
Yo~ ~~ = = nn nr rr nr nr nn nn nee nen eee eee %o 
%Facemaping periodic function for 1-Dimension and used in Thesis 1D 
%Upwind code. 
%Written by Benjamin Davis Created: November 2014 
% Department of Applied Mathematics 
% Naval Postgraduate School 
% Monterey, CA 93943-5216 
Yor == = nnn nn nn rr nn nn nn en eee eee eee eee %o 


function [fmp] = facemap (Ne) 
fmp = zeros(2,Ne); 
for e = 1:Ne 


fmp(l1,e) = e; 
fmp(2,e) = e+1; 
if e==Ne 
fmp(2,e) = fmp(1,1); 
end 


end 
end 





Ti 
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APPENDIX C: 
Two-Dimensional Discontinuous Galerkin 





The following are the main MATLAB codes for the two-dimensional problem from 
Chapter 4. 








%This is the main driver for approximating the 2D Acoustic Wave 
%Equation with an upwind flux on a washer grid. The washer grid 
%can be scaled to various sizes based off of user input. 
%Written by Benjamin Davis and Asst. Professor Jeremy Kozdon 


Yo Department of Applied Mathematics 
%o Naval Postgraduate School 
% Monterey , CA 93943-5216 


*%Synopsis: Discontinuous Galerkin Method for wave equation in second 
*order form using inexact integration and an upwind flux. The outputs 
Yare currently five plots: Convergence rates for w and P and plots of 
Ythe numerical and exact solutions. 


% ee ee ee ee eS ee ee ee ee ee ee ae ee ee ee ee, ee ee ee % 
clear 
N 2% 
n= 1; 
e= 1; 


t_final = pi/2; 
%Define beta to be zero and will run with no theta dependence. 
beta = 4; 


2/% skew_mesh = 0 rectangular mesh 
% skew_mesh = 1 skew element (straight sided) 
% skew_mesh = 2. :: skew element (curved elements) 


skew_mesh = 2; 

%Define anonymous functions for the solution 
rex = @x,y) sqrt(x42 + y%2); 

theta = @(x,y) acos(x/r_ex(x,y)); 

g_ex = @r_ex) besselj(beta,r_ex); 

f_ex = @x,y,t) sin(cs#t-beta»theta(x,y)); 
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P_ex = @(x,y,t) f_ex(x,y,t)*#g_ex(r_ex(x,y)); Y%p(x,y,t) 

w_ex = @x,y,t) cxcos(c+t-betastheta(x,y))*g _ex(r_ex(x,y)); %dp/dt 
%Select inner and outer radius of washer. 

%Found these values using fzero command on Bessel function. 

rl = fzero(g_ex,6); 

r2 = fzero(g_ex,14); 

disp ([rl,r2]); 


%Full Washer. 
Yrn is what 2*pi will be divided by to change the size of the mesh. 
%For example, if you use 4, you will get pi/2 or 1/4 of the washer. 


rm = 4; 


%Scaling the washer. 
Yrn = ceil ((2+pi+r2)/(r2-r1)); 


%For storing information, Don't Change 
z=1; 
for nel = 2.%(1:3) 


nq = N+1; %Number of integration/quadrature points 
ngl = N+1; %Number of interpolation points in one direction 
nelx = nel; 


nely = nel; 


%lnterpolation and Integration Points 

[psx,w] = legendre_gauss_lobatto (N+1); 

%Create Grid 

[coord ,intma ,bsido ,iperiodic ,Np,Ne,nboun,nface] = create_grid_2d ( 
nelx ,nely ,N, psx, plot_grid ,skew_mesh) ; 


fprintf('Number of Elements %4d with polynomial order %2d\n' ,Ne,N) 


%Create Sides/Edge Information for DG 

[iside ,jeside] = create_side(intma,bsido ,Np,Ne,nboun, nface ,ngl) ; 
[face ,imapl,imapr] = create_face(iside ,intma,nface ,ngl) ; 

s_face = face; 

%Changing face code to enforce new BC. 

[face] = faceBC(nface, face); 


80 





74 


78 


84 


88 


90 


92 


94 


96 


98 


100 


102 
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%Will not be periodic with face from faceBC. If you want periodic, 
Ycomment out faceBC function. 

face = create_face_periodicity (iside , face ,coord ,nface ,nboun) ; 
%Building the washer. 

[Rad] = radius(coord ,r1,r2); 
[theta] = thetapolar(coord,rn) ; 

[coord] = newcoord (Rad, theta ,coord) ; 

%Construct Lagrange Basis and Jacobian Matrix 

[L,dL] = lagrange_basis(psx,psx) ; 

[ksi_x ,ksi_y ,eta_x ,eta_y ,x_ksi,x_eta,y_ksi,y_eta,jac] = metrics2 ( 





coord ,intma,L,dL,Ne,ngl ,nq) ; 

%Building and Element to Element function. 
%EtoEfunctBC incorporates boundary conditions. 
EtoEBC = EtoEfunctBC (nface , face ,Ne) ; 

%Constuct M, D, D_ksi, D_eta Matrices 

[M,D] = mass_diff2D(L,dL,w) ; 

M1D =M; 

M = kron(M,M) ; 


D_eta = kron(D,eye(ngl)); 
D_ksi = kron(eye(ngl) ,D); 


%Construct L1,L2,L3,L4 
e0 = sparse(1,1,1,ngl,1); 


en = sparse(ngl,1,1,ngl,1); 
L1 = kron(e0',eye(ngl)); 
L2 = kron(eye(ngl) ,en') ; 
L3 = 3 


kron(en',eye(ngl) 
L4 = kron(eye(ngl) ,e0') ; 

%Building Matrix Terms and new Jacobian for solving the RHS for 
equation 5, 6 and 7. 

[MAT, MAT_ksi, MAT_eta, Je ,A,B,C,D,E,F,G,H] = MatrixTerms2D(Ne, D_ksi, 
D_eta, y_eta ,y_ksi,x_eta,x_ksi,jac ,M); 

%Building Dx and Dy for RHS 

[Dx,Dy] = DxDy(Ne, ksi_x , ksi_y ,eta_x ,eta_y , D_ksi, D_eta) ; 

%Building the surface jacobians for the faces 

[Sj1 ,Sj2,5j3 ,Sj4] = SurfaceJac2D(Ne,L1,L2,L3,L4,x_ksi,x_eta,y_ksi, 











81 





106 


108 





20 





124 


30 
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y_eta); 

YCompute the Normals per element 
[nx_1,nx_2,nx_3,nx_4,ny_1,ny_2,ny_3,ny_4] = ElementNormals2D (Ne,L1, 
L2,L3,L4,x_ksi,x_eta,y_ksi,y_eta ,Sj1 ,Sj2 ,Sj3 ,Sj4); 

%Big Matrix of Ones 

M1 = ones(ngl«ngl) ; 








%Building Initial Condition 
for e = 1:Ne 
for i = 1:ngl 


for j 1:ngl 

I = intma(e,i,j); 
x = coord(I,1); 
= coord(I,2); 
qO{1l}(e,i,j) = w_ex(x,y,0); %dp/dt =w at t=0 
qO{2}(e,i,j) = P_ex(x,y,0); %P(x,y,t) at t=0; 


end 

end 
end 
ql = q0; 
Yo VSLSYSISISLSISIS/S/SISIS/S/SISIS/S/S/SIS/S/S/S1SYSISIS/S/S/SIS/S/S/SIS/S/S/SIS/S/S/SYS/SISIS/S/S/SIS/S/SISIS/S/SISIS/SISISLS/SISI0 
Yo. VSLSYSISISLSLSIS/S/S/SIS/S/S/S/S/S/S/SIS/S/S/SISVS/SIS/S/S/SIS/S/S/SIS/S/S/SIS/S/S/S/S/S/SIS/S/S/S/S/S/S/SIS/S/SISIS/SISIS/S/S/SI0 
% Time Step Calculation 
dt = (1/2) «((1/Ne) /(N‘2)); 
Nt = round(t_final/dt); 
dt = t_final/Nt; 


of SHSYSSYSSS SSS SSS SSSSSSSSSSS SSS SSSS SSS SSS SSS SS SSS SSSSSSSSSSSSSSSSSS 
0 /0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0 


Yo. Y1SISIS1S15/5/8/5/8/81915151515/5/5/819151815151515/5/5/5/81915151818/5/5/831915181515/515/5/8/81315151818/8/8/8/818/8/0/018/8/5/o 
% Beginning of the RK54 Time integrator 
a=[ 0.0, 

-567301805773.0/ 1357537059087.0, 

-2404267990393.0/ 2016746695238.0, 

-3550918686646.0/ 2091501179385.0, 

-1275806237668.0/ 842570457699.0]; 


b = [1432997174477.0/ 9575080441755.0, 
5161836677717.0/ 13612068292357.0, 
1720146321549.0/ 2090206949498.0, 
3134564353537.0/ 4481467310338.0, 
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pa 


46 





2277821191437.0/ 14882151754819.0]; 


R{1} = 0; 
R{2} = 0; 
for k = 1:Nt 
qm = q0; 
for s = 1:5 
R{1} = s)#«R{1}; 
R{2} = s)«R{2}; 


%Building all the flux information with Boundary Conditions 

[fw1,fw2,fw3,fw4,fpl,fp2,fp3,fp4] = UpwindFlux2DFwFpBC (Ne,qm 
, L1,L2,L3,L4,EtoEBC,Dx,Dy,nx_1,ny_1,nx_2,ny_2,nx_3,ny_3,nx_4,ny_4); 

%*Solving Right Hand Side 

[R1] = RHSDG2Dn(Ne,MAT, Je ,M,M_1D,M1,qm,Dx,Dy,L1,L2,L3,L4, 
nx_1,ny_1,nx_2,ny_2,nx_3,ny_3,nx_4,ny_4,5j1 ,Sj2 ,Sj3 ,Sj4 ,fw1,fw2,fw3, 
fw4,fp1,fp2,fp3,fp4,c,MAT_ksi, MAT_eta) ; 








Ww 
N 
Il 
wW 
ja 
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+ 
wW 
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end 


q0 = qm; 
end 


®BVSSYS SYS SSSS SSS SSSYSS SSS SSSS SSS SSSSSSSS SSS SSS SSSSSSS SSS SYS SS SSS SASS SSS 
5 


0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0 
PR S/SSISYSYS/S/S1 51 S1S1S1S13/5/5/5/5 51 5S181513/5/5/5/5/S181518/5/5/5/5/8/S{SS151518/5/5/5/5/S1815/8/5/5/S(5/S/S{S(S18/8/8/5/S1 51510 
%For plotting purposes 
%Solving the exact solution at final time. 
for e = 1:Ne 
for i = 1:ngl 
for j = 1:ngl 
I = intma(e,i,j); 
x = coord(I,1); 
y coord(I,2); 
qex{1}(e,i,j) = w_ex(x,y,t_final); 


qex{2}(e,i,j) = P_ex(x,y,t_final); 
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178 end 

end 

180 end 

%Plotting the exact solution vs. numerical solution. 


182 m=1; 
for e = 1:Ne 
184 for i = 1:ngl 
for j = 1:ngl 

186 

w(m) = qm{1}(e,i,j); 
188 wexact(m) = qex{1}(e,i,j); 
190 pm) = qm{2}(e,i,j); 

pexact(m) = gex{2}(e,i,j); 
193 

I = intma(e,i,j); 
194 xp(m) = coord(I); 

man+1; 
196 end 

end 

198 end 


%Building the Error Plot Information 
200 L2errw = 0; 
L2errp = 0; 


202 


for elm = 1:Ne 


204 pnts = (elm-1)+#(N+1)42 + (1:(N+1)%2); 
dw = w(pnts)-wexact(pnts) ; 
206 dp = p(pnts)-pexact(pnts) ; 
L2errw = dwsMs Je {elm}+dw'+L2errw ; 
208 L2errp = dpsMsJe {elm}+dp'+L2errp ; 
end 


L2err(1,z) = sqrt(Np); 


216 





L2err(2,z) = 
L2err(3,z) = 
z= zt; 


%Convergence 


rate_w = (log(L2err(2,1:end-1))-log(L2err(2,2:end))) 


sqrt (L2errw) ; 
sqrt(L2errp) ; 


Rates 
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./ (log (L2err 





(1,2:end))-log(L2err(1,l:end-1))); 
rate_p = (log(L2err(3,1:end-1))-log(L2err(3,2:end))) ./ (log(L2err 
(1,2:end))-log(L2err(1,l:end-1))); 


fprintf('Convergence rate for w: ') 
220 disp (rate_w) 


222 fprintf('Convergence rate for p: ') 
disp (rate_p) 
224 end 


226 | YEYSYSYS/S/S/5/8/5/815815/5/5/8/81815/8/5/5/5/5151515/5/5/8/5/8151515/8/8/81818/35/5/5/8/813/3/5/5/8/8/8 1515/0 


PSYSSSSY SSS SISSY SSYSS SSS SS SISSY SSS SSS SISSY SSS SSS SV SSSISSYS 
0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0 


228|% il 2 3 4 5 6 7 8 
c = [1 0 0;0 1 0;0 0 0; 0 0 0; 00 0; 00 0; 0 0 0; 0.1 0.2 0.8;... 
230 00 0;0 0 0;0 0 0; 00 0; 00 0; 00 0; 00 0; 0.5 0.3 0.2]; 
% 9 10 11 12 13 14 15 16 


o32| figure (1) ; 

WS/S//S/S/S/S/S/S/S/SISAN Convergence P10 t%s75/18/5/8/5/8/5/8/5/8/5/8/8/0 

24] subplot(1,2,1); %hold on 

loglog(L2err(1,:) ,L2err(2,:),'*-','color',c(N,:) ) 


236 title (‘Convergence Plot for w') 
xlabel('Np') 

238 ylabel('L2 w') 
axis tight 


240| YYSS/SIS/S/S/S/S/SIS/S/. Convergence Plo HYSSSSISISISS/SIS/S/S/5/0 
subplot(1,2,2) ;% hold on 


242 loglog(L2err(1,:) ,L2err(3,:),'#-','color',c(N,:)) 
title ('Convergence Plot for p') 

244 xlabel( 'Np') 
ylabel('L2 p') 

246 axis tight 


%Plot Numerical Solution P(x,y,t) 

as] figure (2) ; 

title ('Numerical Solution for P') 

250| Xlabel( 'x') 

ylabel('y') 

252 Zlabel( ‘t') 

plot3(coord(intma,1), coord(intma,2), qm{2}(:),'#') 
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254/%Plot Exact Solution of P(x,y,t) 
figure (3); 
256 title ('Exact Solution for P') 
plot3(coord(intma,1), coord(intma,2), qex{2}(:),'*') 
258)%Plot Numerical Solution w 
figure (4); 
260 title ('Numerical Solution for w') 
% xlabel('x') 
202/% ylabel('y') 
% zlabel('t') 
21| plot3 (coord(intma,1), coord(intma,2), gqm{1}(:),'*') 
%Plot Exact Solution of w 
266) figure (5) ; 
title ('Exact Solution for w') 


268 


co 





plot3(coord(intma,1), coord(intma,2), qex{1}(:),'«') 











nN 


%This is the main driver for approximating the 2D Acoustic Wave 
%Equation with Upwind Flux on a Square Grid. The grid can be rotated 


» 


Yvarious degrees and skewed based on user input. 
%Written by Benjamin Davis and Asst. Professor Jeremy Kozdon 


6| Yo Department of Applied Mathematics 
%o Naval Postgraduate School 
3] % Monterey, CA 93943-5216 


%poSynopsis: Discontinuous Galerkin Method for wave equation in 





o/%second order form using inexact integration and upwind flux. The 
“outputs are currently four plots: Convergence rates for w and P and 
2;%plots of the numerical and exact solutions. 

Aso cst eesed ete eeeetersew eet acesse ses esate seen b eer teobeerseeees aed %o 
4| clear 

6|N 6; 

n= 1; 

sjc = 1; 

t_final = .25; 
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a 


i) 


26 


34 


38 


40 


44 


48 


54 





% skew_mesh = 0 :: rectangular mesh 


ll 
= 


% skew_mesh skew element (straight sided) 
% skew_mesh = 2. :: skew element (curved elements) 


skew_mesh = 2; 


%Rotation Angle 
grid_rotation_angle=0; “CCW rotation in degrees 
%For information storage. Don't change. 
z=1; 
for nel = 2.%(1:3) 
nq = N+1; %Number of integration/quadrature points 
ngl = N+1; %Number of interpolation points in one direction of an 
element 
nelx = nel; 
nely = nel; 


plot_grid = 0; %l will display the grid, 0 will not display the grid 


%lnterpolation and Integration Points 

[psx ,w] = legendre_gauss_lobatto (N+1); 

%Create Grid 

[coord ,intma,bsido ,iperiodic ,Np,Ne,nboun,nface] = create_grid_2d ( 
nelx ,nely ,N, psx, plot_grid ,skew_mesh) ; 


fprintf('Number of Elements %4d with polynomial order %2d\n' ,Ne,N) 
%Rotate Grid 

[coord_rotated] = rotate_grid_v2 (coord ,intma,Np,Ne,ngl, plot_grid, 
grid_rotation_angle); 

%Store Rotated COORDS 

coord=coord_rotated ; 

%Construct Lagrange Basis and Jacobian Matrix 

[L,dL] = lagrange_basis(psx,psx); 

[ksi_x ,ksi_y ,eta_x ,eta_y ,x_ksi,x_eta,y_ksi,y_eta,jac] = metrics2 ( 





coord ,intma,L,dL,Ne,ngl ,nq) ; 

%Create Sides/Edge Information for DG 

[iside ,jeside] = create_side(intma,bsido ,Np,Ne,nboun, nface ,ngl) ; 
[face ,imapl,imapr] = create_face(iside ,intma,nface ,ngl); 

face = create_face_periodicity (iside , face ,coord ,nface ,nboun) ; 
%Building and Element to Element function 
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64 


66 


68 


74 


80 
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EtoE = EtoEfunct(nface , face ,Ne) ; 
%Constuct M, D, D_ksi, D_eta Matrices 
[M,D] = mass_diff2D(L,dL,w) ; 

M_1D =M; 

M = kron(M,M) ; 


D_eta = kron(D,eye(ngl)); 
D_ksi = kron(eye(ngl) ,D); 
%Construct L1,L2,L3,L4 

e0 = sparse(1,1,1,ngl,1); 
en = sparse(ngl,1,1,ngl,1); 


Ll = kron(e0',eye(ngl)); 
L2 '); 
L3 
L4 = kron(eye(ngl) ,e0') ; 


kron(eye(ngl) ,en 


); 
kron(en',eye(ngl)); 
%Building Matrix Terms and new Jacobian for solving the RHS for 
equation 5, 6 and 7. 

[MAT, MAT_ksi, MAT_eta, Je ,A,B,C,D,E,F,G,H] = MatrixTerms2D (Ne, D_ksi, 
D_eta, y_eta ,y_ksi,x_eta,x_ksi,jac ,M); 

%Building Dx and Dy for RHS 

[Dx,Dy] = DxDy(Ne, ksi_x , ksi_y ,eta_x ,eta_y , D_ksi, D_eta) ; 

%Building the surface jacobians for the faces 

[Sj1 ,Sj2,5j3 ,Sj4] = SurfaceJac2D(Ne,L1,L2,L3,L4,x_ksi,x_eta,y_ksi, 
y_eta); 

YCompute the Normals per element 











[nx_1,nx_2,nx_3,nx_4,ny_1,ny_2,ny_3,ny_4] = ElementNormals2D (Ne,L1, 
L2,L3,L4,x_ksi,x_eta,y_ksi,y_eta ,Sjl ,Sj2 ,5j3 ,Sj4); 

%Big Matrix of Ones 

M1 = ones(ngl«ngl) ; 








%*Building Initial Condition 
%Equation for initial condition is the following: 
Ydp/dt = nxpi*xsin(n«pisx)*sin(ns«pixy) «cos (nspi«t) 
%P(x,y,t) = sin(nspi+x)+sin(nspixy)+sin(nspi+t) 
for e = 1:Ne 
for i = 1:ngl 
for j = 1:ngl 
I 


x 


intma(e,i,j); 
coord(I,1) ; 
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00 


02 


04 





06 


08 





Nv 











= coord(I,2); 
qO{1l}(e,i,j) = sqrt(2)+n+pi+sin(nspisx)*sin(nspisy) ; 
O{2}(e,i,j) = 0; 
end 
end 
end 
ql = 40; 
Yo. ASISISIS/S/S/SISISLS/S/S/S/S/S/SISISIS/S/S/S/S/S1S/S/S/S/S/S/SUSLS/S/S/S/S/S/S1S/S/S/S/S/S1S1S/S/S/S/S/S/SISUS/S/S/S/SISISIS/S/S/o 
Yo. YSISISIS/S/S/SISISIS/S/S/S/S/S/SIS/S/515/5/5/5/815/5/5/S/S/8/S1818/5/5/5/8/8/81 51 S15/5/5/81S1S/S/5/5/S/S/81S1 S/S/S/S/S/S1S15/8/S/o 
% Time Step Calculation 
dt = ((1/Ne)/(N‘%2)); 
Nt = round(t_final/dt); 
dt = t_final/Nt; 


of SHYSSSSSSS SSS SASS SSS SSS SSSSSSSSSSS SSS SSSA SSS S SSS SSS SYSSSSSSSSSS SSS, 
0 /0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0 


Yo. Y/SISIS1S15/5/8/5/8181515151515/5/8/819/5151515151515/5/8/81915181518/5/5/819/518151515/8/5/5/8819151515/8/5/5/8/8/8/0/5/8/8/5/o 
% Beginning of the RK54 Time integrator 

[ 0.0, 

-567301805773.0/ 1357537059087.0, 

-2404267990393.0/ 2016746695238.0, 

-3550918686646.0/ 2091501179385.0, 

-1275806237668.0/ 842570457699.0]; 


a 


ion 
ll 


[1432997174477.0/ 9575080441755.0, 
5161836677717.0/ 13612068292357.0, 
1720146321549.0/ 2090206949498.0, 
3134564353537.0/ 4481467310338.0, 
2277821191437.0/ 14882151754819.0]; 


for s = 1:5 
R(1} = a(s)«R{1}; 
R{2} = a(s)s#R{2}; 
%Building all the flux information 
[fw1,fw2,fw3,fw4,fpl1,fp2,fp3,fp4] = UpwindFlux2DFwFp(Ne,qm, 
L1,L2,L3,L4, EtoE ,Dx,Dy,nx_1,ny_1,nx_2,ny_2,nx_3,ny_3,nx_4,ny_4); 





89 














cad 


) 


i 


46 


48 


60 





%*Solving Right hand side 

[R1] = RHSDG2Dn(Ne,MAT, Je ,M,M_1D,M1,qm,Dx,Dy,L1,L2,L3,L4, 
nx_l,ny_1l, nx_2,ny_2,nx_3,ny_3,nx_4,ny_4,5j1 ,Sj2 ,Sj3 ,Sj4 ,fwl1,fw2,fw3 
,fw4,fpl1,fp2,fp3,fp4,c,MAT_ksi, MAT_eta) ; 





end 


q0 = qm; 
end 


®PBVASYS SYS SSSS SS SS SSS SYSSSSSSSSSSSSSSSYSSSSSSSSSSS SS SSS SSS SSSSSSS S/S 
0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0 


PBSSSYS SYS SSS SSSSS SSS SASSI SSS SS SS SSS SSS SSS S SSS SSS SSS SSS SSSSSSES 
0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0 


2|%For plotting purposes 


%*Solving the exact solution at final time. 
for e = 1:Ne 
for i = 1:ngl 

for j = 1:ngl 
I = intma(e,i,j); 
x = coord(I,1); 
y coord(I,2); 
gex{1}(e,i,j) = sqrt(2)+nspissin(nspi+x)»*sin(n#pixy) cos ( 


sqrt(2)snspi+t_final); 
gqex{2}(e,i,j) = sin(n+pi+x)*sin(nspisy)*sin (sqrt (2) #nspix 
t_final); 
end 
end 
end 
%Plotting the exact solution vs. numerical solution. 
m=1; 
for e = 1:Ne 
for i = 1:ngl 
for j = 1:ngl 


w(m) = qm{1}(e,i,j); 


wexact(m) = qex{1}(e,i,j); 
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p(m) = qm{2}(e,i,j); 
pexact(m) = gex{2}(e,i,j); 


I = intma(e,i,j); 
68 xp(m) = coord(I); 
m=mn+1; 
70 end 
end 





72 end 
%*Building the Error Plot Information 


174 L2errw = 0; 
L2errp = 0; 
for elm = 1:Ne 
78 pnts = (elm-1)+#(N+1)42 + (1:(N+1)%2); 
dw = w(pnts)-wexact(pnts) ; 





end 
L2err(1,z) 


80 dp = p(pnts)-pexact(pnts) ; 
L2errw = dwsMs: Je {elm}sdw'+L2errw ; 
82 L2errp = dp«MsJe{elm}«dp'+L2errp; 


= sqrt(Np); 


L2err(2,z) 
L2err(3,z) 
z+1; 
%Convergence Rates 
rate_Ww (log(L2err(2,1: 
(1,2:end))-log(L2err(1,1: 
rate_p = (log(L2err(3,1: 
(1,2:end))-log(L2err(1,1: 


= sqrt(L2errw 


sqrt (L2errp 


Z = 


3 
); 


end-1))-log(L2err(2,2:end) )) 
end-1))); 
end-1))-log(L2err(3,2:end) )) 
end-1))); 


./ (log (L2err 


./ (log (L2err 


fprintf('Convergence rate for w: ') 
disp (rate_w) 
194 
fprintf('Convergence rate for p: ') 
disp (rate_p) 
end 
P/STSISYS/5151513/515151513/815151818/515151518/5151518/3/518181818/5151818/515151313/8/5151818/515/8/5/o 


PBVASYSSYSSYSYSS SSS SSS SSSSS SSS SSS SSS SSSS SSS SS SSSSSSSSSS 
0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0 


% il 2 3 4 5 6 7 8 





200 
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c = [1 0 0;0 1 0;0 0 0; 00 0; 00 0; 00 0; 0 0 0; 0.1 0.2 0.8;... 
202 00 0;0 0 0;0 0 0; 00 0; 00 0; 00 0; 00 0; 0.5 0.3 0.2]; 
% 9 10 Dat 12 13 14 15 16 

204] figure (1); 

Y/SISISIHS/S/S/I//0N Convergence P10 tVs/6/5/5/5/5/5/8/5/8/815/5/5/5/o 

206] subplot(1,2,1); %hold on 

loglog(L2err(1,:) ,L2err(2,:),'color',c(N;:) ) 


208 title ('Convergence Plot for w') 
xlabel( 'Np') 

210 ylabel('L2 w') 
axis tight 


212| YYS/S/SIS/S/SIS/S/S/S/S/. Convergence Plo YSSISISISISISIS/SIS/S/5/5/0 


subplot(1,2,2) ;% hold on 


214 loglog(L2err(1,:) ,L2err(3,:),'color',c(N;,:) ) 
title ('Convergence Plot for p') 

216 xlabel( 'Np') 
ylabel('L2 p') 

218 axis tight 


%Plot Numerical Solution P(x,y,t) 
220] figure (2) ; 
title ('Numerical Solution for P') 
222| Xlabel( 'x') 
ylabel('y') 
24| Zlabel( ‘'t') 
plot3(coord(intma,1), coord(intma,2), qm{2}(:),'«') 
226|%Plot Exact Solution of P(x,y,t) 
figure (3); 
2s) title ('Exact Solution for P') 
plot3(coord(intma,1), coord(intma,2), qex{2}(:),'«') 
230}%Plot Numerical Solution w 


figure (4); 

232 title ('Numerical Solution for w') 
plot3(coord(intma,1), coord(intma,2), qm{1}(:),'*') 
234(%Plot Exact Solution of w 

figure (5); 

236 title ('Exact Solution for w') 





plot3(coord(intma,1), coord(intma,2), qex{1}(:),'#') 
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%o 


3 % 


%o 


Function for computing fw and fw with an Upwind Flux with p=0 
Boundary Conditions enforced on a washer mesh. 
Written by Benjamin Davis and Asst. Professor Jeremy Kozdon 
Department of Applied Mathematics 
Naval Postgraduate School 
Monterey, CA 93943-5216 





woo en ee ee eee - -- - 2 2 -- -- - - 2 2 2-22 - - - 2-2 ~~ ---- % 
function [fwl,fw2,fw3,fw4,fpl1,fp2,fp3,fp4] = UpwindFlux2DFwFpBC(Ne,ql1,L1 
,L2,L3,L4,EtoEBC,Dx,... 
Dy, nx_1,ny_1,nx_2,ny_2,nx_3,ny_3,nx_4,ny_4) 
¢= 1; 
for k = 1:Ne 


y 


nE1 = EtoEBC(k,1) 
nE2 = EtoEBC(k,2) ; 
nE3 = EtoEBC(k,3) ; 
nE4 = EtoEBC(k,4) 


, 


%Side one of element face 
wml = Li+ql{1}(k,:) '; 
wpl = L3%ql{1}(nE1,:) '; 
if k == nEl1 
wpl = -wml; 
end 
%Side two of element face 
wm2 = L2+ql{1}(k,:) '; 
wp2 = L4sql{1}(nE2,:) '; 
if k == nE2 
wp2 = -wm2; 
end 
%Side three of element face 
wn8 = L3xql{1}(k,:) |; 
wp3 = Ll+ql{1}(nE3,:) '; 
if k == nE3 
wp3 = -wm3; 
end 
%Side four of element face 
wn = L4eql{1}(k,:) |; 
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wp4 = L2sql{1}(nE4,:) '; 
if k == nE4 

wp4 = -wm4; 
end 


%Side one of element face 


pml = nx_1{k}+L1sDx{k}«q1{2}(k, 
ppl = nx_3{nE1}+L3+Dx{nE1}«ql1 {2} 


2}(nE1,:) '; 
if k == nEl 
ppl = -pm1; 
end 


%Side two of element face 


%Side Three of element ie 
pm3 = nx_3{k}+L3%Dx{k}«q ser 
pp3 = ae eee ee a 


2}(nE3,:) '; 


if k == nE3 
pp3 = -pm3; 
end 
%Side Four of element face 


pm4 = nx_4{k}+L4sDx{k}«q1{2}(k, 
pp4 = nx_2{nE4}+L2sDx{nE4}+ ql {2 


2}(nE4,:) '; 
if k == nE4; 
pp4 = -pm4; 
end 


>) '+ny_l1{ 


>) '+ny_3{ 


>) '+ny_4{ 
}(nE4,:) '+ny_2{nE4} em 


k}+L1sDy{k 1{2}(k,:) '; 
(nE1,:) '+ny_3{nE1} eae 


pm2 = nx_2{k}+L2sDx{k}+ql{2}(k,:) '+ny_2{k}+L2s«Dy{k 1{2}(k,:) '; 
pp2 = nx_4{nE2}+L4.Dx{nE2}%q1{2}(nE2,:) '+ny_4{nE2} Sone 
2}(nE2,:) '; 
if k == nE2 
pp2 = -pm2; 
end 


}+L3xDy{k}+q1{2}(k,:) '; 
2}(nE3,:) oe eee eee 


k}+L4sDy{k 1{2}(k,:) '; 


Her aaa e of fpl,fp2,fp3,fp4,fw1,fw2,fw3,fw4 


fpl{k} = (1/2)*(pml - ppl) + 1/(2+c) *(wpl-wml) ; 
Aaa = (1/2) *(pm2 - pp2) + 1/(2+c) *(wp2-wmd) ; 
fp3{k} = (1/2) *(pm3 - pp3) + 1/(2c) «(wp3-wms3) ; 
fp4{k} = (1/2) +*(pm4 - pp4) + 1/(2+c) «(wp4-wm4) ; 
fwl{k} = (1/2) +(wpl - wml) - c/2#(ppl+pm1) ; 
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fw2{k} = (1/2) +(wp2 - wm2) - c/2+#(pp2+pm2) ; 
fw3{k} = (1/2) +(wp3 - wm3) - c/2+#(pp3+pm3) ; 
fw4{k} = (1/2) «(wp4 - wm4) - c/2s(pp4tpm4) ; 
end 
end 
%G- ------------ 2-2-2 2 2 2 eee % 


%Function for computing fw and fp with and Upwind flux with no boundary 


3|%conditions being enforced. Used with periodic boundary conditions. 


%Written by Benjamin Davis and Asst. Professor Jeremy Kozdon 





5|%o Department of Applied Mathematics 
Yo Naval Postgraduate School 
7|% Monterey , CA 93943-5216 
Yor = =~ = ne ne re re ee ee ee ee 2-2-2 ----- % 
o|}function [fwl,fw2,fw3,fw4,fpl,fp2,fp3,fp4] = UpwindFlux2DFwFp(Ne,q1,L1, 
L2,L3,L4,EtoE,Dx,... 
Dy, nx_1,ny_1,nx_2,ny_2,nx_3,ny_3,nx_4,ny_4) 
1 eS 1 
for k = 1:Ne 








nEl = EtoE(k,1); 
nE2 = EtoE(k,2) ; 
nE3 = EtoE(k,3) 
nE4 = EtoE(k,4) ; 

%Side one of element face 
wml = Li+ql{1}(k,:) '; 

wpl = L3%ql{1}(nE1,:) '; 
%Side two of element face 
wm2 = L2+ql{1}(k,:) '; 

wp2 = L4sql{1}(nE2,:) '; 
%Side three of element face 
wn8 = L3xql{1}(k,:) |; 

wp3 = Ll«ql{1}(nE3,:) '; 
%Side four of element face 
wn = L4eql{1}(k,:) '; 

wp4 = L2«ql{1}(nE4,:) '; 


%Side one of element face 


y 
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pps = 
{2} (nE3 ,:) 
%Side 
pm4 = 
pp4 = 
{2}(nE4,:) 


%Calculation of fpl,fp2,fp3,fp4,fwl,fw2,fw3,fw4 
= (1/2)*(pml - ppl) + 1/(2+*c)*(wpl-wml) ; 
) ) + 1/(2*c) «(wp2-wmd) ; 
= (1/2) +(pm3 - pp3) + 1/(2sc 
) ) + 1/(2%c) «(wp4-wm4 


fp1{k 
fp2{k 
fp3{k 

{k 


end 


end 


nx_1{k}+L1sDx{k}+q1{2}(k,:) '+ny_1{k}+L1sDy{k}+q1{2}(k,:) '; 
nx_3{nE1}+L3+Dx{nE1}+ql1{2}(nE1,:) '+ny_3{nE1}+L3sDy{nE1}+q1 


'. 
, 


two of element face 


nx_2{k}+L2sDx{k}*«q1{2}(k,:) '+ny_2{k}+L2sDy{k}+q1{2}(k,:) |; 
nx_4{nE2}+L4sDx{nE2}sq1{2}(nE2,:) '+ny_4{nE2}+L4sDy{nE2}+q1 


'. 
, 


Three of element face 


nx_3{k}+L3+Dx{k}«q1{2}(k,:) '+ny_3{k}+L3«Dy{k}+q1{2}(k,:) |; 
nx_1{nE3}+L1+Dx{nE3}+ ql {2}(nE3,:) '+ny_1{nE3}+L1sDy{nE3}+q1 


Vs 
, 


Four of element face 


nx_4{k}+L4sDx{k}+q1{2}(k,:) '+ny_4{k}+L4sDy{k}+q1{2}(k,:) |; 
nx_2{nE4}+L2sDx{nE4}+q1{2}(nE4,:) '+ny_2{nE4}+L2.Dy{nE4}.+q1 


ee 
, 


= (1/2) +(pm2 - pp2 
= (1/2) +(pm4 - pp4 


= (1/2 
= (1/2 
= (1/2 
= (1/2 


+(wpl - wml 
*(wp2 - wm2 
«(wp3 - wm3 
+(wp4 - wm 


Vw eww 


- c/2«(ppl+pm!) ; 
- c/2«(pp2+pm2 
- c/2s(pp3+pm3 
- c/2+(pp4tpm4 











2}%Central Flux Routine. 


%for the two-dimensional problem. 


%Written by Benjamin Davis 


%o 


5| Yo 


%o 


Department of Applied Mathematics 


Naval Postgraduate School 
Monterey , CA 93943-5216 


Provides the center flux values for fw and fp 











function [fwl,fw2,fw3,fw4,fpl,fp2,fp3,fp4] = CenterFlux2DFwFp(Ne,q1,L1, 


L2,L3,L4,EtoE,Dx,... 
Dy, nx_1,ny_1,nx_2,ny_2,nx_3,ny_3,nx_4,ny_4) 
for k = 1:Ne 
nEl = EtoE(k,1); 
nE2 = EtoE(k,2) ; 
nE3 = EtoE(k,3) 
nE4 = EtoE(k,4); 
%Side one of element face 
wml = Li+ql{1}(k,:) '; 
wpl = L3sql{1}(nE1,:) '; 
%Side two of element face 
wm2 = L2+ql{1}(k,:) '; 
wp2 = L4sql{1}(nE2,:) '; 
%Side three of element face 
wnB = L3+ql{1}(k,:) '; 
wps = Lisql{1}(nE3,:) '; 
%Side four of element face 
wnt = L4eql{1}(k,:) '; 
wp4 = L2«q1{1}(nE4,:) '; 





y 


fw1{k} = (1/2) *(wpl - wml) ; 
fw2{k} = (1/2) *(wp2 - wmd); 
fw3{k} = (1/2) +(wp3 - wmd) ; 
fw4{k} = (1/2) *(wp4 - wm); 


%Calculation of fpl,fp2,fp3,fp4 
%Side one of element face 


pml = nx_1{k}+L1sDx{k}«ql{2}(k,:) '+ny_1{k}+L1lsDy{k}+q1{2}(k,:) ' 


y 


ppl = nx_3{nE1}+L3+Dx{nE1}s ql {2}(nE1,:) '+ny_3{nE1}+L3sDy{nE1}+q1 


[2b MEL 2) "2 
%Side two of element face 


pm2 = nx_2{k}+L2sDx{k}*«ql{2}(k,:) '+ny_2{k}+L2«Dy{k}+q1{2}(k,:) ' 


, 


pp2 = nx_4{nE2}«L4+Dx{nE2}+ql {2}(nE2,:) '+ny_4{nE2}«L4«Dy{nE2}+q1 


{2}(nE2,:) '; 
%Side Three of element face 


pm3 = nx_3{k}+L3sDx{k}«ql{2}(k,:) '+ny_3{k}+L3«Dy{k}+q1{2}(k,:) ' 


y 


pp3 = nx_1{nE3}+L1+Dx{nE3}« ql {2}(nE3,:) '+ny_1{nE3}+L1sDy{nE3}+q1 


{2}(nE3,:) '; 
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%Side Four of element face 


pm4 = nx_4{k}+L4sDx{k}+q1{2}(k,:) '+ny_4{k}+L4sDy{k}+q1{2}(k,:) |; 
pp4 = nx_2{nE4}+L2+Dx{nE4}%q1{2}(nE4,:) '+ny_2{nE4}»+L2sDy{nE4}+q1 








48 fpl{k} = (1/2) +(pml1 - ppl); 
fp2{k} = (1/2) +(pm2 - pp2); 
50 fp3{k} = (1/2) +(pm3 - pp3); 
fp4{k} = (1/2) +(pm4 - pp4); 
2 end 
end 
| Yor = - - - 2-2 oo oe ee oo ee ee 2 o  - - 2 ~~~ -- 
% Function to solve the combined discretized equations from Chapter 4 
3|% for the 2D Acoustic Wave equation. Computes RHS and inverts the 
% required matrices to compute w and p for second order wave equation. 
5/% Written by Benjamin Davis and Asst. Professor Jeremy Kozdon 
% Department of Applied Mathematics 
7|%o Naval Postgraduate School 
% Monterey , CA 93943-5216 
(Cater ee ceeeesecaeebereeweesewesseeeeeoseseseeddeecesteedéoedeessessesee 








function [R] = RHSDG2Dn(Ne,MAT, Je ,M,M_1D,M1, qs ,Dx,Dy,L1,L2,L3,L4,nx_1, 
ny_1,nx_2,ny_2,nx_3,ny_3,nx_4,ny_4,Sjl1 ,Sj2 ,5j3 ,Sj4 ,fw1,fw2,fw3,fw4, 
fpl,fp2,fp3,fp4,c,MAT_ksi, MAT eta) 








bal 


= zeros(size(qs{1})); 


w 
N 
Il 


zeros (size(qs{2})); 


for e = 1:Ne 


R{2}(e,:) = (MAT{e}+MI1'sJefe}«M)«qs{1}(e,:) |; 


R{2}(e,:) = R{2}(e,:)' + ((Dx{e} '+L1'snx_l{e}+ Dy{e}'+L1'+ny_1{e}) 


+M_1DsSjl{e}+fwlf{e}); 


R{2}(e,:) = R{2}(e,:)' + ((Dx{e} '+L2'snx_2{e}+ Dy{e}'+L2'+ny_2{e}) 


+M_1DsSj2{e}«fw2{e}) ; 
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R{2}(e,:) = R{2}(e,:)' + ((Dx{e} '+L3'snx_3{e}+ Dy{e}'+L3'+ny_3{e}) 


+M_1DsSj3{e}+fw3{e}) ; 


R{2}(e,:) = R{2}(e,:)' + ((Dx{e} '+L4'snx_4{e}+ Dy{e}'+L4'sny_4{e}) 


#M 1D. Sj4{e}+fw4{e}) ; 


R{1}(e,:) = c42s(L1'sM1DsSjl{e}+fplf{e}+ L2'sM 1DsSj2{e}sfp2{e}+L3 


'«M_1D+ $j3 {e}«fp3 {e}+L4'«M_1Ds Sj4 {e}«fp4{e}) ; 


R{1}(e,:) = R{1}(e,:)' - cA22(MAT_ksi{e}+MAT_eta{e})»qs{2}(e,:) |; 


Mil{e} = (MAT{e}+M1'« Je{e}M) ; 
[L_1,U_1] = lu(Mll{e}, ‘vector'); 
Ble Ua.) (hal) Rize) 3 
R1 = R1'; 


R{2}(e,:) = R1; 


M1_2{e} = Je{e}«M; 
[L_2,U_2] = lu(M1_2{e}, ‘vector'); 


K2 = U2.\ (L2 1 RT} Ce;:) 3 
R2 = R2'; 
R{1}(e,:) = R2; 


end 


3}end 











%This function computes the LGL grid and elements in 2D. Modified by 


Yamesh . 


5 % 


%Function given to F.X. Giraldo's MA4245 class August 2014. 
%Modified: Jan 2015 
%Written by F.X. Giraldo on 4/2008 
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3|%Jeremy Kozdon and Benjamin Davis for thesis project to produce a skewed 











w 
oS 





%o 
%o 
% 
% 
% 
% 
% 
% 
%o 
% 
Yo 
% 
%o 
%o 
% 
% 


5| % 


Department of Applied Mathematics 


Naval Postgraduate School 
Monterey, CA 93943-5216 
INPUT LIST: nelx and nely are the number of elements in x and y 


nop is the polynomial order 


xgl are the interpolation points on the element. 


coord are the coordinates: x=coord(:,1) and y=coord(: ,2) 


intma is the connectivity list that points to the global 


bsido is the boundary data (used by ISIDE and FACE) 


iperiodic points to another point if periodicity is 


of global points 

of elements 

of boundary edges 

of faces/edges in the grid 


OUTPUT LIST: 
gridpoint number 
applicable 
npoin = number 
nelem = number 
nboun = number 
nface = number 


function [coord,intma,bsido ,iperiodic ,npoin,nelem,nboun,nface] = 


create_grid_2d(nelx ,nely ,nop,xgl,plot_grid , skew_mesh) 


% Define Grid Dimensions 


ngl=nop +1; 


npoin=(nopsnelx + 1)#(nopsnely + 1); 


nelem=nelxsnely ; 


nboun=2snelx + 2snely; 


nface=2snelem + nelx + nely; 


%Initialize Global Arrays 


coord=zeros (npoin ,2) ; 


intma=zeros (nelem,ngl ,ngl) ; 


bsido=zeros (nboun,4) ; 


iperiodic=zeros(npoin,1) ; 


%Initialize Local Arrays 


node=zeros (npoin ,npoin) ; 


%Set some constants 


xmin=-1; 


3] xmax=+1; 


ymin = -1; 


5] ymax=+1; 


dx=(xmax-xmin) /nelx ; 
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dy=(ymax- ymin) /nely ; 
nop=ngl -1; 


oJnx=nelxsnop + 1; 


ny=nely+nop + 1; 
YSSENERATE COORD 


ip =0; 
s| jj =0; 
for k=1:nely 
yO=ymin + real(k-1) «dy; 
if (k == 1) 
11=1; 
else 
11 =2; 
end 
for l=l1:ngl 
jj=jjtl, 
ii=0; 


for i=1:nelx 

x0=xmin + real (i-1)+dx; 

xc = xO + [0,1;0,1]«dx; 

ye = yO + [0,0;1,1]«+dy; 

if (skew_mesh==1) 
xc=xc +(1/8)+sin(pisyc) .*(1-xc) .*(1+xc) ;%x; 
yc=yc+(1/8)+sin(pixxc).*«(1-yc).«(1+yc) ;%y; 

end 

if (i == 1) 
jl =1; 

else 


end 
for j=jl:ngl 
ii=ii + 1; 
ip=ip + 1; 
ax=(xgl(j)+1)/2; 
ay=(xgl(1)+1)/2; 
x = xc(1,1)*(1-ax)s(1l-ay) + xc(1,2)+ax«(l-ay)... 
+ xc(2,1)+#(1-ax)+#( ay) + xc(2,2)#axs#( ay); 
y = ye(1,1)+#(1-ax)*#(1-ay) + ye(1,2)saxs(l-ay)... 
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+ yco(2,1)+*(1-ax)+( ay) + ye(2,2)#ax#( ay); 


7 coord (ip ,1)=x; 
coord (ip ,2)=y; 
fe if (skew_mesh==2) 


coord(ip,1)=x +(1/8)*sin(pisy) «(1-x)#(1+x) ;%x; 
91 coord(ip ,2)=y +(1/8)*#sin(pi+x)«(1-y)«(1+y) ;%y; 
end 
node(ii , jj )=ip; 
end %j 


end %i 
end %l 
end %k 
99] “GENERATE INTMA 


ie =0; 


N 


ii} for k=1:nely 

for i=1:nelx 

103 ie=ie+1; 

for l=1:ngl 

105 jj =(mgl-1) #(k-1) + 1; 
for j=1:ngl 

107 ii=(ngl-1)+(i-1) + j; 
ip=node(ii ,jj); 


109 intma(ie,j,1l)=ip; 
end %j 
111 end %l 
end %i 
us}end %k 
%Generate BSIDO 
us| ib =0; 
for i=1:nelx 
117 ie=i; 
ib=ib+1; 
119 il=(i-1)+(ngl-1) + 1; 
i2=(i-1)+(ngl-1) + nel; 
121 ipl=node(il ,1); 
ip2=node(i2 ,1); 
123 bsido(ib ,1)=ip1; 





bsido (ib ,2)=ip2; 
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bsido(ib ,3)=ie; 
bsido (ib ,4) =6; 


27} end 


%Right Boundary 


o|}for i=1:nely 


ie=(nelx) (i); 
ib=ib+1; 
il=(i-1)+(ngl-1) + 1; 
i2=(i-1)+(ngl-1) + nel; 
ipl=node(nx, il); 
ip2=node(nx,i2); 
bsido(ib ,1)=ip1; 
bsido(ib ,2)=ip2; 
bsido(ib ,3)=ie; 
bsido(ib ,4) =6; 

end 

%lop Boundary 

for i=nelx:-1:1 
ie=nelem - (nelx - i); 
ib=ib+1; 
il=(i-1)+(ngl-1) + nel; 
i2=(i-1)+(ngl-1) + 1; 
ipl=node(il ,ny) ; 
ip2=node(i2 ,ny) ; 
bsido(ib ,1)=ip1; 
bsido(ib ,2)=ip2; 
bsido (ib ,3)=ie; 
bsido(ib ,4) =6; 


53| end 


%Left Boundary 


s|for i=nely:-1:1 


ie=(nelx)»(i-1) + 1; 
ib=ib+1; 
il=(i-1)+(ngl-1) + nel; 
i2=(i-1)+(ngl-1) + 1; 
ipl=node(1,il1); 
ip2=node(1,i2); 
bsido(ib ,1)=ip1; 

bsido (ib ,2)=ip2; 
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bsido(ib ,3)= 


ie; 


bsido (ib ,4) =6; 


end 


7|%Periodicity 


for i=1:npoin 


iperiodic(i)= 


end 

%X- Periodicity 

for i=l:ny 
il=node(1,i); 
i2=node(nx,i) 
iperiodic(i2) 

end 


7/%Y- Periodicity 


for i=1:nx 
il=node(i,1) ; 
i2=node(i,ny) 
iperiodic(i2) 
end 


3/%Plot Grid 


i; 


, 


=il; 


y 


=iperiodic(il); 


if (plot_grid == 1) 


x=zeros (5,1) 
y=zeros (5,1) 
figure; 
hold on; 


, 


, 


for e=l1:nelem 


for j=l: 
for 


ngl-1 
i=1:ngl-1 
il=intma(e,i,j); 
i2=intma(e,i+1,j); 
i3=intma(e,i+1,j+1); 
i4=intma(e,i,j+1); 
x(1)=coord (il ,1 ) 
x(2)=coord (i2 ,1 ) 
x(3)=coord(i3,1); y(3)=coord (i3 ,2 
) 
) 


Il 


x(4)=coord(i4,1); y(4)=coord(i4 ,2 
x(5)=coord(i1l,1); y(5)=coord (il ,2 
plot_handle=plot(x,y,'-r'); 

set(plot_handle , 'LineWidth' ,1.5); 
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; y(1)=coord (il ,2); 
; y(2)=coord (i2 ,2); 


) 
) 
) 
) 
) 


y 


, 


y 
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end 

end 
il=intma(e,1,1); 
i2=intma(e,ngl,1) ; 
i3=intma(e,ngl ,ngl); 
i4=intma(e,1,ngl); 
x(1)=coord(il,1); y(1)=coord(il ,2); 
x(2)=coord(i2,1); y(2)=coord(i2 ,2) ; 
x(3)=coord(i3 ,1); y(3)=coord (i3 ,2) ; 
x(4)=coord(i4,1); y(4)=coord(i4 ,2) ; 
x(5)=coord(il ,1); y(5)=coord (il ,2) ; 
plot_handle=plot(x,y,'-b'); 
set(plot_handle , 'LineWidth' ,2); 

end 


title_text=['Grid Plot For: Ne = ' num2str(nelem) ', N = ‘| num2str( 


nop) ]; 
title ([ title_text],'FontSize' ,18); 
xlabel('X','FontSize' ,18); 
ylabel('Y','FontSize',18); 
axis image 
end 











%Function provided to Professor F.X. Giraldo's MA4245 class. 
%Used by Benjamin Davis for his thesis project. 

%This subroutine creates the array ISIDE which stores all of 
Ythe information concerning the sides of all the elements. 
%Written by Francis X. Giraldo on 1/01 


% Naval Postgraduate School 

% Department of Applied Mathematics 

% Monterey , CA 93943-5502 

% INPUT LIST: intma = element connectivity 

Yo bsido = boundary info (which points are on a boundary, 
Yo which element it belongs to and the boundary 
% condition). 

% npoin = number of global points 
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% nelem = number of elements 


Yo nboun = number of boundary faces/edges 

% nside=nface are the number of sides/face/edges in the grid 
Yo ngl = number of points in one direction in an element 

% OUTPUT LIST: iside = face information such as which points are on a 

% face which elements they belong to and, if a 

% boundary, what is the boundary condition. 

% jeside = for each element and each edge gives the FACE 

%o number 

Yo~ - ~~ nn nn nn nn rn nn nn nn nn nn renee ere ee % 
function [iside ,jeside] = create_side(intma,bsido ,npoin ,nelem,nboun, 


nface ,ngl) 
%global arrays 
iside = zeros(nface ,4) ; 
jeside= zeros(nelem,4) ; 
%local arrays 
Iwher = zeros(npoin,1) ; 
Ihowm = zeros(npoin,1) ; 


2)icone = zeros(5snpoin,1) ; 


inode = zeros(4,1); 


s}jnode = zeros(4,1); 


%Fix Inode 
inode (1) =1; 
inode (2)=ngl; 
inode (3)=ngl; 
inode (4) =1; 
jnode (1) =1; 
jnode (2) =1; 
jnode (3)=ngl; 
jnode (4)=ngl; 
%count how many elements own each node 
for in=1:4 
for ie=1:nelem 
ip=intma(ie ,inode(in) ,jnode(in)); 
lhowm(ip)=lhowm(ip) + 1; 
end %ie 
end %in 


%track elements owning each node 


2| lwher (1) =0; 
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for ip=2:npoin 
lwher (ip )=lwher(ip-1) + lhowm(ip -1) ; 
end %ip 
Yanother tracker array 
Ihowm = zeros(npoin,1) ; 
for in=1:4 
for ie=1:nelem 
ip=intma(ie ,inode(in) ,jnode(in)); 
lhowm(ip )=lhowm(ip) + 1; 
jloca=lwher(ip) + lhowm(ip) ; 
icone (jloca)=ie; 
end “%ie 
end %in 


56| Y_OOP OVER THE NODES 


iloca=0; 
for ip=l1:npoin 
ilocl=iloca; 
iele=lhowm(ip) ; 
if (iele ~= 0 ) 
iwher=lwher (ip) ; 
%dOOP OVER THOSE ELEMENTS SURROUNDING NODE IP 
ipl=ip ; 
for iel=l:iele 
ie=icone (iwher+iel ) ; 
%find out position of ip in intma 
for in=1:4 
inl=in; 
ipt=intma(ie ,inode(in) ,jnode(in)); 
if (ipt == ip) 
break 
end 
end %in 
%Check Edge of Element IE which claims IP 
j=0; 
for jnod=1:2:3 
iold=0; 
j=j+l, 
in2=in + jnod; 
if (in2 > 4) 
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92 in2=in2 -4; 


end 
94 ip2=intma(ie ,inode(in2) ,jnode(in2)) ; 
if (ip2 >= ipl) 
96 %check whether side is old or new 
if (iloca ~= ilocl) 
98 for is=ilocl+1l:iloca 


iside(is ,2); 
100 jloca=is; 
if (iside(is ,2) == ip2) 





102 iold=1; 
break; 
04 end 
end %is 
06 end %iloca 
if (iold == 0) 
08 %NEW SIDE 
iloca=iloca + 1; 
10 iside (iloca ,1)=ip1; 
iside(iloca ,2)=ip2; 
112 iside (iloca ,2+j)=ie; 
elseif (iold == 1) 
114 YOLD SIDE 
iside (jloca ,2+j )=ie; 
116 end %iold 
end %ip2 
118 end %jnod 
end %iel 
120 %Perform some Shifting to order the nodes of a side in QW 
direction 


for is=ilocl+l:iloca 


122 if (iside(is ,3) == 0) 
iside(is ,3)=iside (is ,4) ; 
124 iside(is ,4) =0; 
iside(is ,1)=iside(is ,2) ; 
126 iside(is ,2)=ip1; 
end %iside 
128 end %is 





end %if iele 
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end %ip 


if (iloca ~= nface) 
disp( ‘Error in SIDE. iloca nface = '); 
iloca 
nface 
pause 


36] end 


%RESET THE BOUNDARY MARKERS 
for is=1:nface 
if (iside(is ,4) == 0) 
il=iside(is ,1); 
ir=iside (is ,2) ; 
ie=iside(is ,3); 
for ib=1:nboun 
ibe=bsido(ib ,3); 
ibc=bsido (ib ,4) ; 
if (ibe == ie) 
ilb=bsido(ib ,1); 
irb=bsido (ib ,2) ; 
if (ilb == il & irb == ir) 
iside (is ,4)=-ibc; 
break 
end %ilb 
end %ibe 
end %ib 
end %iside 


s6}end %is 


Y%FORM. ELEMENT/SIDE CONNECTIVITY ARRAY 
for is=l:nface 


y 


iel=iside (is ,3 


isl=iside(is ,1); 


y 


) 
ier=iside (is ,4) ; 

) 

) 


is2=iside (is ,2 
% LEFT SIDE 
for in=1:4 


y 


il=intma(iel ,inode(in) ,jnode(in)); 
inl=in + 1; 
if (inl > 4) 


inl=1; 
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end %in1 
i2=intma(iel ,inode(in1) ,jnode(in1)); 
if ((isl1 == il) && (is2 == i2)) 
jeside(iel ,in)=is; 
end %is1 
end %in 
%RIGHT SIDE 
if (ier > 0) 
for in=1:4 
il=intma(ier ,inode(in) ,jnode(in)) ; 
inl=in + 1; 
if (inl > 4) 
inl=1; 
end %in1 
i2=intma(ier ,inode(inl) ,jnode(in1)); 
if ((isl == i2) && (is2 == il1)) 
jeside (ier ,in)=is; 











end %is1l 
end %in 
end %ier 
end %is 
Vee e See se Se SSeS SSSSe SeSSees SSeS Se See Skee SSSeee Se SeSeus Sees eee esse Se % 


%Function provided to Professor F.X. Giraldo's MA4245 class. Used by 
%Benjamin Davis. This subroutine constructs the Side Information for 
%a High Order Spectal Element Quads 

%Written by Francis X. Giraldo 


%o Department of Applied Mathematics 

% Naval Postgraduate School 

% Monterey , CA 93943-5216 

% INPUT LIST: iside = face information to know which points are on a 
% face and which elements they belong to 

% intma = element connectivity 

Yo nface = number of faces/edges in the grid 

Yo ngl = number of interpolation points in one direction 
% in an element 
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N 


3) face=zeros(nface ,4) ; 


5}imapr=zeros (4,2,ng 


5}% OUIPUT LIST: face = face information stating which local edge number 


Yo the face is on and which elements they belong to 

Yo (left and right neighbors).are the metric terms 

% needed to imapl and imapr give the tensor-product (i,j) 
9|% points of the edge points. 
a % 


function [face ,imapl,imapr]=create_face(iside ,intma,nface ,ngl) 
%global arrays 


imapl=zeros (4,2,ngl 
1); 


3 
) 


%local arrays 


7|inode=zeros (4,1); 


jnode=zeros (4,1); 


9|%Construct ihe Pointer 


inode(1)= 
si inode (2)= a 
inode (3)=ngl; 
33) inode (4) =1; 
jnode (1) =1; 
35] jnode (2) =1; 
jnode (3)=ngl; 


N 


w 


41 





jnode (4)=ngl; 
%Construct IMAP arrays 


o}for l=1:ngl 


%eta=-1 
imapl(1,1,1)=1; 
imapl(1,2,1)=1; 
imapr(1,1,1)=ngl+1-1; 
ke ,2,1)=1; 


Yksi=+ 

or 1,1)=ngl; 
imapl(2,2,1)=1; 
imapr(2,1,1)=ngl; 
imapr(2,2,1)=ngl+1-1; 


%eta=t+1 
imapl(3,1,1)=ngl+1-1; 
imapl(3,2,1)=ngl; 
imapr(3,1,1)=1; 
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imapr(3,2,1)=ngl; 


%ksi=-1 
imapl(4,1,1)=1; 
imapl(4,2,1)=ngl+1-1; 
imapr(4,1,1)=1; 
imapr(4,2,1)=1; 


end %l 
%loop thru the sides 
for i=1:nface 


y 


ipl=iside(i,1 


iel=iside(i,3); 


y 


) 
ip2=iside(i,2); 
) 
ier=iside(i,4); 
%check for position on Left Element 
for j=1:4 
jl=j; 
j2=j+1; 
if (j2 > 4) 
j2=1; 
end %j2 
jpl=intma(iel ,inode(j1),jnode(j1)); 
jp2=intma(iel ,inode(j2),jnode(j2)); 
if (ipl == jpl && ip2 == jp2) 
face(i,1)=j; 
break; %leave J loop 
end %ip1 
end %j 
%check for position on Right Element 
if (ier > 0) 
for j=1:4 
ji=j; 
j2=j+1; 
if (j2 > 4) 
j2=1; 
end %j2 
jpl=intma(ier ,inode(jl) ,jnode(j1)); 
jp2=intma(ier ,inode(j2) ,jnode(j2)); 


if (ipl == jp2 && ip2 == jp1) 
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93 face (i,2)=j; 


break; leave J loop 
95 end %ip1 
end %j 
97 end %ier 


%Store Elements into face 
99 face(i,3)=iel; 
face(i,4)=ier; 

iifend %i 














%Function for changing the face code to enforce p=0 boundary conditions 
ging Pp y 


w 


%for 2nd order acoustic wave equation on a washer. 
%Written by Benjamin Davis Created: 09 May 2015 


5|% Department of Applied Mathematics 

Yo Naval Postgraduate School 

71% Monterey , CA 93943-5216 

% ake peels agua geet at ee Caran eee em oe age eS ah a, eee meen eed a oe See ea ee rl ge ae! eh ae bees ot ay TE ee a ee aol et ee SOE TiS Oe ee eye ee, PY RT ee % 


oO 


function[face] = faceBC(nface , face) 
for k = 1:nface 


1 if face(k,1) == 2 
if face(k,4) == -6 
3 face(k,4) = -1; 
end 
5 end 
end 
7 for k = 1:nface 
if face(k,1) == 
9 if face(k,4) == -6 





face(k,4) = -1; 
21 end 

end 

23 end 

end 
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%Function provided by Professor F.X. Giraldo to his MA4245 class. 
%Used by Benjamin Davis 
%This subroutine builds Periodic BCs along the 4 edges of a rectangular 


%domain . 

%Written by Francis X. Giraldo on 2/2007 

Yo Naval Postgraduate School 

%o Department of Applied Mathematics 

% Monterey, CA 93943-5216 

% INPUT LIST: iside = face information 

%o face = more face information 

% coord = gridpoint coordinates 

% nface = number of faces 

% nboun = number of boundaries 

% OUTPUT LIST: face = augments the FACE data structure to include 

%o periodicity 

Yo ~~ = nn rn nn nr nn rn nn nn nn nnn eer eee Yo 
function face = create_face_periodicity (iside , face ,coord ,nface ,nboun) 
%Constant 

tol=le-6; 


%Local arrays 
ileft=zeros (nboun/4 ,1) ; 
iright=zeros (nboun/4 ,1) ; 
itop=zeros (nboun/4 ,1) ; 
ibot=zeros (nboun/4 ,1); 
%initialize 
nleft=0; nright=0; ntop=0; nbot=0; 
%Find Extrema of Domain 
xmax=max(coord(:,1)); xmin=min(coord(:,1)); 
ymax=max(coord(:,2)); ymin=min(coord(: ,2)); 
%loop thru sides and extract Left, Right, Bot, and Top 
for is=1:nface 
%Check for Periodicity Edges 
ier=face(is ,4); 
if (ier == -6) 
il=iside(is,1); i2=iside(is ,2); 
xm=0.5%( coord(il ,1) + coord(i2,1) ); 
ym=0.5( coord(il ,2) + coord(i2,2) ); 
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%check Grid Point 
if ( abs(xm - xmin) < tol ) %left boundary 
nleft=nleft + 1; 
ileft(nleft)=is; 
elseif ( abs(xm - xmax) < tol ) %right boundary 
nright=nright + 1; 
iright(nright )=is; 
elseif ( abs(ym - ymin) < tol ) %bottom boundary 
nbot=nbot + 1; 
ibot (nbot)=is ; 
elseif ( abs(ym - ymax) < tol ) %top boundary 
ntop=ntop + 1; 
itop (ntop)=is ; 
else 
disp('No match in PERIODIC_BCS for is ier = '); 
is 
ier 
pause 
end %if 
end %ier 
end %is 
%Loop through Periodic BCs 
%First: Do Left and Right 


2}for i=1:nleft 


isl=ileft(i); 
il=iside(isl ,1); 
yll=coord (il ,2); 
%*Search for Corresponding Right Edge 
for j=1:nright 
isr=iright(j); 
i2=iside(isr ,2); 
yr2=coord (i2 ,2) ; 
if ( abs(yll-yr2) < tol ) %they match 
face(isl ,2)=face(isr ,1); 
face(isl ,4)=face(isr ,3) ; 
face(isr ,3)=-6; Y%means skip it due to Periodicity 
iside(isl ,4)=iside(isr ,3); 
iside(isr ,3) =-6; 
break; 
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78 end %if 

end %j 

sojend %i 

%Second: Do Top and Bottom 

s2)for i=1:nbot 

isl=ibot(i); 

84 il=iside(isl ,1); 

xll=coord(il ,1); 

86 %Search for Corresponding Top Edge 
for j=1:ntop 


88 isr=itop(j); 
i2=iside(isr ,2); 
90 xr2=coord(i2 ,1) ; 
if ( abs(xll-xr2) < tol ) %they match 
92 face(isl ,2)=face(isr ,1); 
face(isl ,4)=face(isr ,3); 
94 face(isr ,3) =-6; Y%means skip it due to Periodicity 
iside(isl ,4)=iside(isr ,3); 
96 iside(isr ,3) =-6; 
break; 
98 end %if 
end %j 


1oo}end %i 














Nv 


%Function given to Professor F.X. Giraldo's MA4245 class. Modified by 
%Ben Davis and Jeremy Kozdon. 


r 


%This function computes the Metric Terms for General 2D Quad Grids. 
%Written by F.X. Giraldo on 4/2008 





6|% Department of Applied Mathematics 
Yo Naval Postgraduate School 
8} % Monterey , CA 93943-5216 
% INPUT LIST: coord = gridpoint coordinates 
10| % intma = element connectivity 
% psi = basis functions defined as psi(NGL,NQ) 
12|% dpsi = derivative of basis functions defined as 
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% dpsi (NGL,NQ) 


%o nelem = number of elements 

Yo ngl = number of interpolation points in one direction in 
% an element 

% nq = number of integration/quadrature points in one 

% direction in an element. 


% OUTPUT LIST: ksi_x ,ksi_y ,eta_x ,eta_y 








are the metric terms needed to 


% compute derivatives in 

% physical space 

% x_ksi,x_eta,y_ksi,y_eta = are the metric terms needed to 
Yo compute derivatives in 

% physical space 

Yo jac = weightssJacobian defined as jac (nelem,nq,nq) 
A %o 


function [ksi_x,ksi_y ,eta_x,eta_y ,x_ksi,x_eta,y_ksi,y_eta,jac] = 





metrics2 (coord ,intma, psi ,dpsi ,nelem ,ngl ,nq) 


%Initialize Global Arrays 
ksi_x=zeros (nelem ,nq,nq) ; 
ksi_y=zeros (nelem ,nq,nq) ; 
eta_x=zeros (nelem,nq,nq) ; 
eta_y=zeros (nelem,nq,nq) ; 
jac=zeros(nelem,nq,nq) ; 
%I Initialize Local Arrays 
x_ksi=zeros (nelem,nqg,nq) ; 
x_eta=zeros (nelem,nq,nq) ; 
y_ksi=zeros (nelem ,nq,nq) ; 
y_eta=zeros (nelem ,nq,nq) ; 
x=zeros(ngl,ngl); 
y=zeros(ngl,ngl); 
%loop thru the elements 
for ie=1:nelem 
%Store Element Variables 
for j=1:ngl 
for i=1:ngl 
ip=intma(ie,i,j); 
x(i,j)=coord(ip,1); 
y(i,j )=coord (ip ,2) ; 
end %i 
end %j 
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%Construct Mapping Derivatives: dx/dksi, dx/deta, dy/dksi, dy/deta 
[x_ksi(ie ,: ,:) ,x_eta(ie ,:,:)]=map_deriv(psi,dpsi,x,ngl,nq) ; 
[y_ksi(ie,:,:) ,y_eta(ie ,:,:) ]=map_deriv(psi,dpsi,y,ngl,nq) ; 
%Construct Inverse Mapping: dksi/dx, dksi/dy, deta/dx, deta/dy 


for j=1:ng 

for i=1:ng 
xjac=x_ksi(ie,i,j)*y_eta(ie,i,j) - x_eta(ie,i,j)*y_ksi(ie,i,j); 
ksi_x (ie ,i,j)=+1.0/xjac+y_eta(ie,i,j); 
ksi_y(ie,i,j)=-1.0/xjac+x_eta(ie,i,j); 
eta_x(ie,i,j)=-1.0/xjac+y_ksi(ie,i,j); 
eta_y(ie,i,j)=+1.0/xjac+x_ksi(ie,i,j); 


jac(ie,i,j)=abs(xjac); 
end %i 
end %j 
end %ie 











% Function created for taking the x coordinates from coord and 


3/% turning them into radius polar coordinates based on user inputs. 


%Written by Benjamin Davis %Created: 14 April 2015 


5|% Department of Applied Mathematics 


% Naval Postgraduate School 
% Monterey, CA 93943-5216 


function[R] = radius(coord,rl1,r2) 
n = length(coord(:,1)); 
for i= 1:n 
R(i,1) = ((r2-1r1)/2)+coord(i,1) + (rl + 12)/2; 
end 
end 
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%Function to turn the y coord into theta. Conversion from Cartesian to 
polar coordinates. Also, rn is a scaling portion of the washer. 
%Written by Benjamin Davis Created: 14 April 2015 


Yo Department of Applied Mathematics 

% Naval Postgraduate School 

% Monterey , CA 93943-5216 

WG--- - 2 oe en nn ee eee eee eee eee ee --------- % 
function [theta] = thetapolar(coord ,rn) 


n = length(coord (: ,2)); 


t1 =- 0; 
r2 = (2spi)/rn; 
for i= I:n 


theta(i,1) = ((r2-r1)/2)*coord(i,2) + (rl + r2)/2; 
end 
end 











%Function to overwrite original coord for the conversion to polar mesh. 
%Function is required for the washer mesh. 
%Written by Benjamin Davis Created: 14 April 2015 


% Department of Applied Mathematics 

%o Naval Postgraduate School 

% Monterey , CA 93943-5216 

% a ear eae et es ed tem ee eat See ea et epee ar ae ey eee ee ee ne Rae See ne et pee cree tae a ae) ed ee ee me ae Pam ee ee ean eae Mee eae Se ep ee er ery ee eg ee ees ee aed oe ae eee es, ee neg eS ee % 


function[coord] = newcoord(R, theta ,coord) 
n = length(coord(:,1)); 
for i= 1:n 
coord(i,1) = R(i)«cos(theta(i)); 
coord (i ,2) R(i)+#sin(theta(i)); 


end 
end 
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%Element to Element function with Boundary Conditions for Washer Grid 
*Synopsis: Element to element function for identifying and storing 
%information for looking across the faces of an element to it's 
‘corresponding neighbor. 

%Written by Benjamin Davis 


Yo Department of Applied Mathematics 

%o Naval Postgraduate School 

% Monterey , CA 93943-5216 

% 

Yo pL = Face of Left Element 

Yo pR = Face of Right Element 

Yo Ls = Number of Left Element 

% Rs = Number of Right Element 

%Output : EtoE is a Matrix of Number of Elements (Row) vs. the Number 


%of Sides per element (Four for a square grid). For each element, it 
%displays its corresponding neighbor with respect to boundary 
%conditions . 


function [EtoEBC] = EtoEfunctBC (nface , face ,Ne) 
EtoEBC = zeros (Ne,4) ; 
for l=1:nface 
pL = face(1,1); %Face of left element 


pR = face(1,2); %Face of Right Element 
Ls = face(1,3); %Number of Left Element 
Rs = face(1,4); %Number of Right Element 
if Rs == -1 
EtoEBC(Ls,pL) = Ls; 
elseif (Ls ~= -6 && Rs ~= -6) 


EtoEBC(Ls,pL) = Rs; 
EtoEBC(Rs,pR) = Ls; 
end 
end 


5}end 
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%Element to Element function for a square grid with periodic boundary 
*conditions . 

YoSynopsis: Element to element function for identifying and storing 
%information for looking across the faces of an element to it's 


w 


a 


%corresponding neighbor. 


N 


%Written by Benjamin Davis 


%o Department of Applied Mathematics 

9|% Naval Postgraduate School 

Yo Monterey, CA 93943-5216 

1|% 

Yo pL = Face of Left Element 

3|% pR = Face of Right Element 

Yo Ls = Number of Left Element 

5|% Rs = Number of Right Element 

Yo 

7|%Output : EtoE is a Matrix of of Number of Elements (Row) vs. the 
Number 


%*of Sides per element (Four for a square grid). For each element, it 
%displays its corresponding neighbor. 


oO 


2}function [EtoE] = EtoEfunct(nface , face ,Ne) 
EtoE = zeros(Ne,4) ; 

23 for l=1:nface 

pL = face(1,1); 


25 pR = face(1,2); 
Ls = face(1,3); 

27 Rs = face(1,4); 

29 if (Ls ~= -6 && Rs ~= -6) 

EtoE(Ls,pL) = Rs; 

31 EtoE(Rs,pR) = Ls; 
end 

33 end 

end 
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%Function for building the Mass and Differentiation matrices for 
%2-Dimension . 
%Written by Benjamin Davis Created: July 2014 Modified Jan 2015 


Yo Department of Applied Mathematics 

% Naval Postgraduate School 

% Monterey, CA 93943-5216 

% a a a a % 


function [M,D] = mass_diff2D(L,dL,w) 
n = length(L(:,1)); 

M = zeros(n,n); 

D = zeros(n,n); 











for i = lin 
for j = 1:n 
M(i,j) = ((LG,:) #LQ,:) )ew(1,:)"); 
DGi,j) = LCi ,:)*dLQj,:)'; 
end 
end 
end 
A %o 
%Written by Benjamin Davis 
Yo Department of Applied Mathematics 
% Naval Postgraduate School 
% Monterey, CA 93943-5216 
% 23: Jan 2015 


*Synopsis: Building Matrix Terms in order to solve the equations for 
%the 2D Acoustic Wave Equation. 
%Output: Matrix elements used to solve the RHS in the RK54 scheme. 


function [MAT,MAT_ksi, MAT_eta, Je ,A,B,C,D,E,F,G,H] = MatrixTerms2D (Ne, 
D_ksi, D_eta, y_eta, y_ksi ,x_eta ,x_ksi ,jac ,M) 
for k = 1:Ne 
Je{k} = diag(jac(k,:)); 
Jeinv{k} = diag (1./jac(k,:)); 
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A = D_ksi's(diag(y_eta(k,:) )«Jeinv {k}«Msdiag(y_eta(k,:) ))+D_ksi; 
B = D_eta'»(diag(y_ksi(k,:) )«Jeinv {k}»Msdiag (y_eta(k,:) ))+D_ksi; 
C = D_ksi'+(diag(y_eta(k,:) )*Jeinv {k}«M:+ diag (y_ksi(k,:) ))*D_eta; 
D = D_eta'+(diag(y_ksi(k,:) )#Jeinv {k}sMsdiag (y_ksi(k,:) ))+D_eta; 
E = D_ksi'+(diag(x_eta(k,:) )sJeinv {k}sMsdiag(x_eta(k,:) ))*D_ksi; 
F = D_eta'+(diag(x_ksi(k,:) )*Jeinv {k}*Ms diag (x_eta(k,:) ))*D_ksi; 
G = D_ksi'+(diag(x_eta(k,:) )#Jeinv {k}sMsdiag(x_ksi(k,:) ))+D_eta; 
H = D_eta's(diag(x_ksi(k,:) )+Jeinv {k}«Ms diag (x_ksi(k,:) ))+D_eta; 
MAT{k} = A-B-C+D+E-F-G+H; 
MAT _ksi{k} = A-B+E-F; 
MAT _eta{k} = -C+D-GtH;%%%C+D-G+H 
end 
end 
Go- ~~ oe ne re rn rn oo ee ee ee ee ~~ -- % 


%Function to build Dx and Dy used to solve the RHS in the RK54 scheme. 
%Written by Benjamin Davis 





% Department of Applied Mathematics 
% Naval Postgraduate School 
51% Monterey, CA 93943-5216 
%o 23 Jan 2015 
Yor == = nnn nr nn rn nn nn nn nnn nee eee eee eee %o 
function [Dx,Dy] = DxDy(Ne, ksi_x , ksi_y , eta_x ,eta_y , D_ksi, D_eta) 
for k = 1:Ne 

Dx{k} = diag(ksi_x(k,:))+*D_ksi + diag(eta_x(k,:) )+D_eta; 


= 
a 
I 


diag(ksi_y(k,:) )*D_ksi + diag(eta_y(k,:) )+D_eta; 


end 
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2;%Function for building the Surface Jacobians for the faces on a 2D grid 


% Written by Benjamin Davis 


%o Department of Applied Mathematics 

%o Naval Postgraduate School 

% Monterey , CA 93943-5216 

%o 23 Jan 2015 

Yor = = = ww nn nn en nn nn en en ne ee eee =e % 


function [Sjl ,Sj2,$j3,Sj4] = SurfaceJac2D(Ne,L1,L2,L3,L4,x_ksi,x_eta, 
y_ksi, y_eta) 
for k = 1:Ne 
%Building the surface jacobians for the faces 
Sj1_3{k} = sqrt(x_ksi(k,:) ..2 + y_ksi(k,:) .%2); 
Sj2_4{k} = sqrt(x_eta(k,:) .42 + y_eta(k,:) .%2); 
%Isolate Surface Jacobians individually for four faces 
Sjl{k} = diag(L1+Sj1_3{k}'); 














Sj3{k} = diag(L3+Sj1_3{k}'); 
Sj2{k} = diag(L2+Sj2_4{k}'); 
Sj4{k} = diag(L4+Sj2_4{k}'); 
end 
end 
% ee ee ee ee ee een ee ee ee ee a a a a ee, ae ae eee ae ee ee ee ey ae a ee eee eke a a eee pe oe ae ey, oe ee ee % 


%Function to build the element normals for the fours faces on a 2D grid 
%Written by Benjamin Davis 


% Department of Applied Mathematics 

%o Naval Postgraduate School 

% Monterey, CA 93943-5216 

%o 23 Jan 2015 

Yor - = 2 wn nn nn nn nn nn nn en rn nn ne ne ne nee eee ee % 


function [nx_1,nx_2,nx_3,nx_4,ny_1,ny_2,ny_3,ny_4] = ElementNormals2D (Ne 
,L1,L2,L3,L4,x_ksi,x_eta ,y_ksi,y_eta,Sj1 ,Sj2 ,Sj3 , Sj4) 
for k=1:Ne 


yCompute the Normals per element 
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nx_l{k} = diag(Ll+y_ksi(k,:) './diag(Sj1l{k})); 
ny_l{k} = -diag(L1+x_ksi(k,:) './diag(Sj1{k})); 


nx_2{k} = diag(L2sy_eta(k,:) './diag(Sj2{k})); 
ny_2{k} = -diag(L2sx_eta(k,:) './diag(Sj2{k})); 


nx_3{k} = -diag(L3+y_ksi(k,:) './diag(Sj3{k})); 
ny_3{k} = diag(L34x_ksi(k,:) './diag(Sj3{k})); 


nx_4{k} = -diag(L4+y_eta(k,:) './diag(Sj4{k})); 
ny_4{k} = diag(L4+x_eta(k,:) './diag(Sj4{k})); 











end 
end 
(rs ee seer daaeseeee Saeed nee sRee Skee esa eee Rene are eKeneee ees aeeseeeser es Yo 
%Function provided to Professor F.X. Giraldo's MA4245 class. 
%Modified by Asst. Professor Jeremy Kozdon and Benjamin Davis 
%This function rotates the grid and plots it. 
%Written by F.X. Giraldo on 4/2008 
%o Department of Applied Mathematics 
Yo Naval Postgraduate School 
% Monterey , CA 93943-5216 
% 
% INPUT LIST: coord are the coordinates 
% intma is the connectivity list 
%o npoin are the number of global points 
Yo nelem are the number of elements 
% ngl is the number of interpolation points in an element 
Yo (polynomial order + 1)\ 
Yo plot_grid is a switch to either plot or not plot 
%o grid_rotation_angle is the grid rotation in degrees 
% OUTPUT LIST: 
% coord are the new rotated coordinates: x=coord(:,1) 
% and y=coord (: ,2) 
Yor so scr c ccc crc ee ccc n ee secs nee scence esc e esse cess cece s seer nse sscsces Yo 
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plot_grid , grid_rotation_angle) 


nop=ngl -1; 


coord_rotated=zeros(npoin ,2) ; 


%Rotate Grid 


alpha=grid_rotation_anglespi/180; 


for i=1:npoin 


xn=cos(alpha)scoord(i,1) - sin(alpha)scoord (i ,2) ; 


yn=sin(alpha)scoord(i,1) + cos(alpha)»coord(i,2) ; 


coord_rotated (i ,1)=xn; 


coord_rotated (i ,2)=yn; 


end 
%Plot Grid 


if (plot_grid == 1) 
x=zeros(5,1); 


y=zeros(5,1); 


figure; 
hold on; 


for e=l1:nelem 


for j=1:ngl-1 


for 


end 
end 


i=1:ngl-1 
il=intma(e,i,j); 
i2=intma(e,i+l1,j); 
i3=intma(e,i+1,j+1); 
i4=intma(e,i,j+1); 
x(1)=coord_rotated (il ,1 
x(2)=coord_rotated (i2 ,1 
x(3)=coord_rotated (i3 ,1 
x(4)=coord_rotated (i4 ,1 


Ve ew wo 


) 
) 
; y(3) 
) 
) 


plot_handle=plot(x,y,'-r'); 
set(plot_handle , 'LineWidth' ,1.5); 


il=intma(e,1,1); 


i2=intma(e,ngl,1); 


i3=intma(e,ngl,ngl) ; 


i4=intma(e,1,ngl) ; 





x(1)=coord(il,1); y(1)=coord (il ,2) ; 


126 


; y(1)=coord_rotated (il ,2) ; 
; y(2)=coord_rotated (i2 ,2) ; 
coord_rotated (i3 ,2) ; 
; y(4)=coord_rotated (i4 ,2) 
x(5)=coord_rotated(il,1); y(5)=coord_rotated (il ,2) 


2>|}function [coord_rotated] = rotate_grid_v2(coord ,intma,npoin,nelem,ngl, 


y 


y 

















60 x(2)=coord(i2,1); y(2)=coord(i2 ,2) ; 
x(3)=coord(i3,1); y(3)=coord (i3 ,2) ; 
62 x(4)=coord(i4,1); y(4)=coord (i4 ,2) ; 
x(5)=coord(il,1); y(5)=coord(il ,2) ; 
64 plot_handle=plot(x,y,'-b'); 
set(plot_handle , 'LineWidth' ,2); 
66 end 
title_text=['Rotated Grid Plot For: Ne = ' num2str(nelem) ', N = ' 
num2str (nop) ]; 
68 title ([ title_text],'FontSize' ,18); 
70 xlabel('X','FontSize' ,18); 
ylabel('Y','FontSize',18); 
72 axis image 
end 
1] %- ------------------------------- 2-2-2 -- 2-2-2 - 2 -------------- % 


%Function for building the Lagrange Polynomials. 
3/%Written by Benjamin Davis in MA4245 %Created: July 2014 


Yo Department of Applied Mathematics 
5|% Naval Postgraduate School 
Yo Monterey, CA 93943-5216 
7 % See ar ae eh et ete emer en ae ra om aed ae set per ree eae ed gee et et eee te Pee Pome et et eh ae am ae Se Sh he eae eye ne se erm ae oe % % 


function [L,dL] = lagrange_basis(x,z) 


%Nth order interpolation 
n = length(x); 


%Length of the equally spaced grid for k = 1:50 
h = length(z); 

3/%Initialize the Lagrange Matrix 

L = ones(n,h); 

s}dL = zeros(n,h); 

%Computation for Lagrange Matrix 








7 for k = 1:h 
for i= 1:n 
9 for j = 1:n 
dl = 1; 
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21 if j ~= i% If j does not equal i 
%Equation for the Lagrange Polynomial 
23 L(i,k) = (z(k)-x(j)) -/(xG)-xQ)) * LG,k); 
for 1 = lin 
25 if (1 ~= i) & (1 ~= j) 
dl = di«(z(k)-x(1)) ./(x(i)-x(1)); 
27 end 
end 
29 dL(i,k) = dL(i,k) + dl/(x(i)-x(j)); 
end 
31 end 
end 

end 

end 














2}%7VCode given to F.X. Giraldo MA4245 Class July 2014 
%Used by Ben Davis. 
%This code computes the Legendre-Gauss-Lobatto points and weights 


> 


Ywhich are the roots of the Lobatto Polynomials. 
%Written by F.X. Giraldo on 4/2000 
Yo Department of Applied Mathematics 


a 


8} % Naval Postgraduate School 
% Monterey, CA 93943-5216 


function [xgl,wgl] = legendre_gauss_lobatto(P) 
2|p=P -1; 
ph=floor( (p+t1)/2 ); 
for i=1:ph 
x=cos( (2+#i-1)*pi/(2sp+l1) ); 
for k=1:20 
[LO,LO_1,L0_2]=legendre_poly (p,x) ; 
8 %Get new Newton Iteration 
dx =-(1-x%2)+L0_1/(-2+x+LO_1 + (1-x4%2)+L0_2); 
20 xX=x+dx ; 
if (abs(dx) < 1.0e-20) 


i 


a 
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22 


24 


w 
np 


w 
oo 


N 





N 





break 
end 

end 

xgl(p+2-i)=x; 

wegl(p+2-i) =2/(p+(p+1)+*L0%2) ; 
end 
%Check for Zero Root 
if (p+1 ~= 2«ph) 

x=0; 

[LO,LO_1,L0_2]=legendre_poly(p,x); 

xgl(ph+1)=x; 

wegl(ph+1) =2/(p+(p+1)+L0%2) ; 
end 


%Find remainder of roots via symmetry 


6|)for i=1:ph 


xgl(i)=-xgl(p+2-i); 
wel(i)=+wel(pt+2-i); 
end 











%Code given to F.X. Giraldo MA4245 Class July 2014 


3|%Used by Ben Davis 


%This code computes the Legendre Polynomials and its 1st and 2nd 


5}%derivatives 


% 

%Written by F.X. Giraldo on 4/2000 

%o Department of Applied Mathematics 

%o Naval Postgraduate School 

% Monterey, CA 93943-5216 
RR Yo 


function [LO,L0_1,L0_2] = legendre_poly(p,x) 


3) L1=0;L1_1=0;L1_2=0; 


LO=1;L0_1=0;L0_2=0; 


s}for i=l:p 


L2=L1;L2_1=L1_1;L2_2=L1_2; 
L1=L0;L1_1=L0_1;L1_2=L0_2; 
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3}end 


a=(2#i-1)/i; 
b=(i-1)/i; 
LO=axx+Ll1 - b+L2; 


LO_l=a»(L1 + x*L1_1) - bxL2_1; 


LO_2=a+(2+L1_1 + x+L1_2) 


- bx L2_2; 
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